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Preface 


This volume on Advances in Energy System Optimization contains a selection 
of peer-reviewed papers related to the 2nd International Symposium on Energy 
System Optimization (ISESO 2018). The symposium was held at Karlsruhe Institute 
of Technology (KIT) under the symposium theme “Bridging the Gap Between 
Mathematical Modelling and Policy Support” on October 10-11, 2018. ISESO 
2018 was organized by KIT (Institute for Industrial Production (IIP) and Institute 
of Electric Energy Systems and High-Voltage Technology (IEH)), the Heidelberg 
Institute for Theoretical Studies (HITS), Heidelberg University (Engineering Math- 
ematics and Computing Lab (EMCL)), German Aerospace Center (DLR, Institute 
of Engineering Thermodynamics, Department of Energy Systems Analysis), and 
the University of Stuttgart (Institute for Building Energetics, Thermotechnology, 
and Energy Storage). 

The organizing institutes are engaged in a number of collaborative research 
activities aimed at combining interdisciplinary perspectives from mathematics, 
operational research, energy economics, and electrical engineering to develop new 
approaches to integrated energy system and grid modeling designed to solve real- 
world energy problems efficiently. The symposium was now held for the second 
time. By design, ISESO is limited in size as to ensure a productive atmosphere 
for discussion and reflection on how to tackle the many challenges facing today’s 
and tomorrow’s energy systems. Around 50 international participants attended 17 
international presentations from both industry and academia including 2 keynote 
presentations and 15 contributed papers in 6 sessions. The sessions focused on 
diverse challenges in energy systems, ranging from operational to investment 
planning problems, from market economics to technical and environmental con- 
siderations, from distribution grids to transmission grids, and from theoretical 
considerations to data provision concerns and applied case studies. The presenta- 
tions were complemented by a panel discussion on the symposium theme “Bridging 
the Gap Between Mathematical Modelling and Policy Support” involving senior 
experts from academia and industry. 


vi Preface 


The papers in this volume are broadly structured according to the sessions within 
the symposium as outlined below: 


e Optimal Power Flow 

e Energy System Integration 

e Demand Response 

e Planning and Operation of Distribution Grids. 


The editors of this volume served as the organizing committee. We wish to thank all 
the reviewers as well as all the individuals and institutions who worked hard, often 
invisibly, for their tremendous support. In particular, we wish to thank Nico Meyer- 
Hübner and Nils Schween for the coordination of the local organization. Finally, we 
also wish to thank all the participants, speakers, and panelists for their contributions 
to making ISESO a success. 
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Part I 
Optimal Power Flow 


Feasibility vs. Optimality in Distributed A 
AC OPF: A Case Study Considering eed 
ADMM and ALADIN 


Alexander Engelmann and Timm Faulwasser 


Abstract This paper investigates the role of feasible initial guesses and large con- 
sensus-violation penalization in distributed optimization for Optimal Power Flow 
(OPF) problems. Specifically, we discuss the behavior of the Alternating Direction 
of Multipliers Method (ADMM). We show that in case of large consensus-violation 
penalization ADMM might exhibit slow progress. We support this observation by 
an analysis of the algorithmic properties of ADMM. Furthermore, we illustrate our 
findings considering the IEEE 57 bus system and we draw upon a comparison 
of ADMM and the Augmented Lagrangian Alternating Direction Inexact Newton 
(ALADIN) method. 


Keywords Distributed optimal power flow - ALADIN - ADMM 


1 Introduction 


Distributed optimization algorithms for AC Optimal Power Flow (OPF) recently 
gained significant interest as these problems are inherently non-convex and 
often large-scale; i.e. comprising up to several thousand buses [1]. Distributed 
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optimization is considered to be helpful as it allows splitting large OPF problems into 
several smaller subproblems; thus reducing complexity and avoiding the exchange 
of full grid models between subsystems. We refer to [2] for a recent overview of 
distributed optimization and control approaches in power systems. 

One frequently discussed convex distributed optimization method is the Alternat- 
ing Direction of Multipliers Method (ADMM) [3], which is also applied in context 
of AC OPF [1, 4, 5]. ADMM often yields promising results even for large-scale 
power systems [4]. However, ADMM sometimes requires a specific partitioning 
technique and/or a feasible initialization in combination with high consensus- 
violation penalization parameters to converge [1]. The requirement of feasible 
initialization seems quite limiting as it requires solving a centralized inequality- 
constrained power flow problem requiring full topology and load information 
leading to approximately the same complexity as full OPF. 

In previous works [6-8] we suggested applying the Augmented Lagrangian 
Alternating Direction Inexact Newton (ALADIN) method to stochastic and deter- 
ministic OPF problems ranging from 5 to 300 buses. In case a certain line-search 
is applied [9], ALADIN provides global convergence guarantees to local minimizers 
for non-convex problems without the need of feasible initialization. The results in 
[6] underpin that ALADIN is able to outperform ADMM in many cases. This comes 
at cost of a higher per-step information exchange compared with ADMM and a more 
complicated coordination step, cf. [6]. 

In this paper we investigate the interplay of feasible initialization with high 
penalization for consensus violation in ADMM for distributed AC OPF. We illustrate 
our findings on the IEEE 57-bus system. Furthermore, we compare the convergence 
behavior of ADMM to ALADIN not suffering from the practical need for feasible 
initialization [9]. Finally, we provide theoretical results supporting our numerical 
observations for the convergence behavior of ADMM. 

The paper is organized as follows: In Sect. 2 we briefly recap ADMM and ALADIN 
including their convergence properties. Section 3 shows numerical results for the 
IEEE 57-bus system focusing on the influence of the penalization parameter p on 
the convergence behavior of ADMM. Section 4 presents an analysis of the interplay 
between high penalization and a feasible initialization. 


2 ALADIN and ADMM 


For distributed optimization, OPF problems are often formulated in affinely coupled 
separable form 


min Fi) subjectto x; e Ai, Vie R and Yo Ax =0, (I) 
ee ieR ieR 

where the decision vector is divided into sub-vectors x! = its Be on e R™, 

R is the index set of subsystems usually representing geographical areas of a power 

system and local nonlinear constraint sets X; := {x; € R™i | h,(&;) < 0}. 
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Algorithm 1 ADMM 
0 


Initialization: Initial guesses z;, no for alli € R, parameter p. 
Repeat (until convergence): 


1. Parallelizable Step: Solve for alli € R locally 


k aemin weap © ae ra 
x; = argmin Si Qi) + (Aj) Aixi + 7 Aji zi) 3 st. hj) < 0. (2) 
Xi 
2. Update dual variables: 
a = vk + pA; (xk = zy, (3) 
3. Consensus Step: Solve the coordination problem 
min D £ Ax; Aj Ai Axi + AT A; Axi (4a) 
iER 
s.t. > Aa + Ax;) =0. (4b) 
ieR 


4. Update variables: ‘+! — x* + Axt. 


Throughout this work we assume that f; and h; are twice continuously differentiable 
and that all X; are compact. Note that the objective functions f; : Ri — R 
and nonlinear inequality constraints h; : R”! — R”' only depend on x; and that 
coupling between them takes place in the affine consensus constraint ) eR Aixi = 
0 only. There are several ways of formulating OPF problems in form of (1) differing 
in the coupling variables and the type of the power flow equations (polar or 
rectangular), cf. [4, 6, 10]. 

Here, we are interested in solving Problem (1) via ADMM and ALADIN summa- 
rized in Algorithms 1 and 2 respectively.!’* At first glance it is apparent that ADMM 
and ALADIN share several features. For example, in Step (1) of both algorithms, 
local augmented Lagrangians subject to local nonlinear inequality constraints A; 
are minimized in parallel.” Observe that while ADMM maintains multiple local 
Lagrange multipliers A;, ALADIN considers one global Lagrange multiplier vector 
A. In Step (2), ALADIN computes sensitivities B;, g; and C; (which often can 
directly be obtained from the local numerical solver without additional computation) 
whereas ADMM updates the multiplier vectors Àj. 


'We remark that there exist a variety of variants of ADMM, cf. [3, 11]. Here, we choose the 
formulation from [9] in order to highlight similarities between ADMM and ALADIN. 

Note that, due to space limitations, we describe the full-step variant of ALADIN here. To obtain 
convergence guarantees from a remote starting point, a globalization strategy is necessary, cf. [9]. 
3For notational simplicity, we only consider nonlinear inequality constraints here. Nonlinear 
equality constraints g; can be incorporated via a reformulation in terms of two inequality 
constraints, i.e. 0 < gj(x;) < 0. 
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Algorithm 2 ALADIN (full-step variant) 


Initialization: Initial guess (2°, +9), parameters X; > 0, p, u. 
Repeat (until convergence): 


1. Parallelizable Step: Solve for alli € R locally 


2 
“a-l, st. Aj(xi) <0 |é. 
i 


xf = argmin fi) + A)" Aixi + £ 
Xi 


2. Compute sensitivities: Compute Hessian approximations BK, gradients gk and Jacobians of the 
active constraints C es cf. [9]. 
3. Consensus Step: Solve the coordination problem 


: . eT T 
min [24x] Bt Axi + gi Ax} +a! s+&1sl3 


Ax;,s 4 
iER 

s.t. > Ai (xt +Axj)=s [a® 
iER (5) 
Ct Ax; =0 VIER. 


4. Update variables: zkt! — xk + Axt, At! < Q. 
Similarity to ADMM: Remove C* in (5), set BF = pAj Aj, g} = A} AK and 5; = A} Aj, set 
u = œ (i.e. s = 0) and use dual ascent step (3) for updating A* in (4), cf. [9]. 


In Step (3), both algorithms communicate certain information to a central entity 
which then solves a (usually centralized) coordination quadratic program. However, 
ALADIN and ADMM differ in the amount of exchanged information: Whereas 
ADMM only communicates the local primal and dual variables x; and àj, ALADIN 
additionally communicates sensitivities. This is a considerable amount of extra 
information compared with ADMM. However, there exist methods to reduce the 
amount of exchanged information and bounds on the information exchange are 
given in [6]. Another important difference is the computational complexity of the 
coordination step. In many cases, the coordination step in ADMM can be reduced 
to an averaging step based on neighborhood communication only [3], whereas 
in ALADIN the coordination step involves the centralized solution of an equality 
constrained quadratic program. 

In the last step, ADMM updates the primal variables z;, while ALADIN addi- 
tionally updates the dual variables A. Differences of ALADIN and ADMM also 
show up in the convergence speed and their theoretical convergence guarantees: 
Whereas ALADIN guarantees global convergence and quadratic local convergence 
for non-convex problems if a certain globalization strategy is applied [9], few results 
exist for ADMM in the non-convex setting. Recent works [12, 13] investigate the 
convergence of ADMM for special classes of non-convex problems; however, to the 
best of our knowledge OPF problems do not belong to these classes. 
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3 Numerical Results 


Next, we investigate the behavior of ADMM for large p and a feasible initialization to 
illustrate performance differences between ADMM and ALADIN. We consider power 
flow equations in polar form with coupling in active/reactive power and voltage 
angle and magnitude at the boundary between two neighbored regions [6]. We 
consider the IEEE 57-bus system with data from MATPOWER and partitioning as 
in [6] as numerical test case. 

Figures 1 and 2 show the convergence behavior of ADMM (with and without 
feasible initialization (f.)) and ALADIN for several convergence indicators and two 
different penalty parameters p = 10+ and p = 106.4 Therein, the left-handed 
plot depicts the consensus gap || Ax* ||, representing the maximum mismatch of 
coupling variables (active/reactive power and voltage magnitude/angle) at borders 
between two neighbored regions. The second plot shows the objective function value 
f K and the third plot presents the distance to the minimizer Ix* — x*||oo over the 
iteration index k. The right-handed figure shows the nonlinear constraint violation 
Ig(z* )lloo after the consensus steps (4) and (5) of Algorithms 1 and 2 respectively 
representing the maximum violation of the power flow equations.” 

In case of small » = 10*, both ADMM variants behave similar and converge 
slowly towards the optimal solution with slow decrease in consensus violation, 
nonlinear constraint violation and objective function value. If we increase p to 
p = 10° with results shown in Fig. 2, the consensus violation I Ax*||oo gets smaller 
in ADMM with feasible initialization. The reason is that a large p forces ey being 
close to z leading to small || Ax* || as we have Az* = 0 from the consensus step. 
But, at the same time, this also leads to a slower progress in optimality f* compared 
to p = 10+, cf. the second plot in Figs. | and 2. 

On the other hand, this statement does not hold for ADMM with infeasible 
initialization (blue lines in Figs.1 and 2) as the constraints in the local step 


10° Ñ 10° WA — | 
x 

8 = 8 1 
= *s i old 
3 | : 2 10°74 
= ' 2 |[—-ADMM T i = ` 
i iamm tj) ' i 
} - - ALADIN 10° ‘ 
O 0 i 1019 — 

0 50 100 0 50 100 0 50 100 0 50 100 


Fig. 1 Convergence behavior of ADMM with infeasible initialization, ADMM with feasible initial- 
ization (f.) for o = 10* and ALADIN 


4The scaling matrices X; are diagonal. They are chosen to improve convergence. Hence, entries 
corresponding to voltages and phase angles are 100, entries corresponding to powers are set to 1. 
>The power flow equations for the IEEE 57-bus systems are considered as nonlinear equality 
constraints g;(x;) = 0. Hence, g;(x;) # 0 represents a violation of the power flow equations. 
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10* 
10° 5% 0 En] 
imwana els 10° 

Mare eee ki : SSS 
P i 
' R 1 | 
2 N 3 z i 28 1 
2 5 | 4 2 h 8 i a a 
8 10 N | | 1 ® 405} 
= “i 2 = H = N 
= \ | |—ADMM © ' = N 
' 4| [ADM (f)} | ' i 
i -- ALADIN 10°) ' 
1010 — 0 1010 — 

0 50 100 0 50 100 0 50 100 0 50 100 


Fig. 2 Convergence behavior of ADMM with infeasible initialization, ADMM with feasible initial- 
ization (f.) for p = 10° and ALADIN 


= — 10° 
rf 
i 
8 ! 8 
‘h H N 05 
x j 1 = 10 
Pe acy. 4 | |ADMM (f) 
! -- ALADIN 10°} 5 
1010 — 0 fi 4010 
0 50 100 0 50 100 0 50 100 0 50 100 


Fig. 3 Convergence behavior of ADMM with infeasible initialization, ADMM with feasible initial- 
ization (f.) for p = 10!? and ALADIN 


and the consensus step of ADMM enforce an alternating projection between the 
consensus constraint and the local nonlinear constraints. The progress in the 
nonlinear constraint violation ||g(z*)||oo supports this statement. In its extreme, this 
behavior can be observed when using p = 10!* depicted in Fig. 3. There, ADMM 
with feasible initialization produces very small consensus violations and nonlinear 
constraint violations at cost of almost no progress in terms of optimality. 

Here the crucial observation is that in case of feasible initialization and large 
penalization parameter 0 ADMM produces almost feasible iterates at cost of slow 
progress in the objective function values. From this, one is tempted to conclude that 
also for infeasible initializations ADMM will likely converge, cf. Figs. 1, 2, and 3. 
This conclusion is supported by the rather small 57-bus test system. However, it 
deserves to be noted that this conclusion is in general not valid, cf. [1]. 

For ALADIN, we use p = 10° and u = 10’. Comparing the results of 
ADMM with ALADIN, ALADIN shows superior quadratic convergence also in case 
of infeasible initialization. This is inline with the known convergence properties 
of ALADIN [9]. However, the fast convergence comes at the cost of increased 
communication overhead per step, cf. [6] for a more detailed discussion. Note that 
ALADIN involves a more complex coordination step, which is not straightforward 

to solve via neighborhood communication. Furthermore, tuning of p and u can be 
difficult. 

In power systems, usually feasibility is preferred over optimality to ensure a 
stable system operation. Hence, ADMM with feasible initialization and large pọ can 
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in principle be used for OPF as in this case violation of power flow equations and 
generator bounds are expected to be small and one could at least expect certain 
progress towards an optimal solution. Following this reasoning several papers 
consider || A(x* — z*)||, < € as termination criterion [1, 4]. However, as shown in 
the example above, if p is large enough, and ADMM is initialized at a feasible point, 
this termination criterion can always be satisfied in just one iteration if p is chosen 
sufficiently large. Consequently, it is unclear how to ensure a certain degree of 
optimality. An additional question with respect to ADMM is how to obtain a feasible 
initialization. To compute such an initial guess, one has to solve a constrained 
nonlinear least squares problem solving the power flow equations subject to box 
constraints. This is itself a problem of almost the same complexity as the full OPF 
problem. Hence one would again require a computationally powerful central entity 
with full topology and parameter information. Arguably this jeopardizes the initial 
motivation for using distributed optimization methods. 


4 Analysis of ADMM with Feasible Initialization for  — oo 


The results above indicate that large penalization parameters p in combination 
with feasible initialization might lead to pre-mature convergence to a suboptimal 
solution. Arguably, this behavior of ADMM might be straight-forward to see from an 
optimization perspective. However, to the best of our knowledge a mathematically 
rigorous analysis, which is very relevant for OPF, is not available in the literature. 


Proposition 1 (Feasibility and o > oo imply oe _ x e null(A;)) 
Consider the application of ADMM (Algorithm 1) to Problem (1). Suppose that, for 
all k € N, the local problems (2) have unique regular minimizers xk. 6 Fork € N, 


let ak be bounded and, for alli € R, z € Xi. Then, the ADMM iterates satisfy 
lim x*(p) —x* € null(Aj), Vk > &. 
poo 


Proof The proof is divided into four steps. Steps 1-3 establish technical properties 
used to derive the above assertion in Step 4. 
Step 1. At iteration k the local steps of ADMM are 


: . AT p |" 
x (o) = argmin fi(x;) + (at) Aixi + 3 la (x = z) (6) 
xe% 


2 


Now, by assumption all fis are twice continuously differentiable (hence a 
on A; ), 2% is bounded and all z% € Xi. Thus, for alli € R, un n xh (p) = = z + vr 


with ur € null(A;). 


6 A minimizer is regular if the gradients of the active constraints are linear independent [14]. 
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Step 2. The first-order stationarity condition of (6) can be written as 


-VAAD = yf Vh aP = ATA} + pA] Ai (xf - 2), m 


where yT is the multiplier associated to h;. Multiplying the multiplier update 
formula (3) with AJ from the left we obtain AT AKT! = ATA‘ + pAT Aj (xk — 
z’). Combined with (7) this yields ATA! = -vf(ak) — yk vaj(x8). By 
differentiability of f; and hi, compactness of A; and regularity of xk this implies 
boundedness of Ane 


Step 3. Next, we show by contradiction that Axt € null(A;) for all i € R and 
p > oo. Recall the coordination step (4b) in ADMM given by 


‘ pP k k 
min 2 Zar Ai Ai Axi + AE A; Axi s.t. > A + Ax;) = 0. (8) 
Le 1E 


Observe that any Axk e null(A;) is a feasible point to (8) as Vier Aixt = 0. 
Consider a feasible candidate solution Ax; ¢ null(A;) for which Vier Axt + 
Ax;) = 0. Clearly, a A; Axi (p) will be bounded. Hence for a sufficiently large 
value of p, the objective of (8) will be positive. However, for any Ax; € null(A;) 
the objective of (8) is zero, which contradicts optimality of the candidate solution 
Axi ¢ null(A;). Hence, choosing p sufficiently large ensures that any minimizer 
of (8) lies in null(A;). 


Step 4. It remains to show a 


i 


xh + Axk. Given Steps (1)-(3) this yields zk+l = zk + vr + Axk and hence 


22 2 42 
| A; (x — a) = la (x; = z$) 


2 2 2 


= a, In the last step of ADMM we have z*+! = 


u fa (xi -zf ++ Ark) 


Observe that this implies that, for pọ — oo, problem (6) does not change from step 
k to k + 1. This proves the assertion. a 


Corollary 1 (Deterministic code, feasibility, » —> oo implies xt! = x£) 
Assuming that the local subproblems in ADMM are solved deterministically; 


i.e. same problem data yields the same solution. Then under the conditions of 


si š i k 
Proposition | and for p — œ, once ADMM generates a feasible point x; to 


Problem (1), or whenever it is initialized at a feasible point, it will stay at this point 
for all subsequent k > k. 


The above corollary explains the behavior of ADMM for large p in combination 
with feasible initialization often used in power systems [1, 4]. Despite feasible 
iterates are desirable from a power systems point of view, the findings above imply 
that high values of p limit progress in terms of minimizing the objective. 
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Remark I (Behavior of ALADIN for p — œœ) Note that for 0 — oo, ALADIN 
behaves different than ADMM. While the local problems in ALADIN behave similar 
to ADMM, the coordination step in ALADIN is equivalent to a sequential quadratic 
programming step. This helps avoiding premature convergence and it ensures 
decrease of f in the coordination step [9]. 


5 Conclusions 


This method-oriented work investigated the interplay of penalization of consensus 
violation and feasible initialization in ADMM. We found that—despite often working 
reasonably with a good choice of p and infeasible initialization—in case of feasible 
initialization combined with large values of 0 ADMM typically stays feasible yet 
it may stall at a suboptimal solution. We provided analytical results supporting 
this observation. However, computing a feasible initialization is itself a problem 
of almost the same complexity as the full OPF problem; in some sense partially 
jeopardizing the advantages of distributed optimization methods. Thus distributed 
methods providing rigorous convergence guarantees while allowing for infeasible 
initialization are of interest. One such alternative method is ALADIN [9] exhibiting 
convergence properties at cost of an enlarged communication overhead and a more 
complex coordination step [6]. 
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Security Analysis of Embedded HVDC A 
in Transmission Grids PEE 


Marco Giuntoli and Susanne Schmitt 


Abstract Security-constrained optimal power flow is widely used by grid operators 
during their short-term operations. Possible contingencies (e.g. outages of lines, 
generators) analyzed during this phase may be very critical for the stability of 
the power system. Thus, a simplified or reduced approach might not be desired. 
Moreover, the presence of embedded high-voltage direct current (HVDC) links (i.e. 
the DC link is done between two or more synchronous AC nodes) increases the 
degrees of freedom of the physical problem and thus an economically better working 
point for the system can be reached. To do this, a coordinated control between 
the AC and DC systems is required using one single stage of optimal power flow 
including the security assessment. In this work the SC-OPF (Security-Constrained 
Optimal Power Flow) approach is used to analyse the effects of contingency cases 
within the DC system. 


Keywords Optimal power flow - Security-constrained optimal power flow - 
Embedded HVDC - N-1 criterion 


1 Introduction 


There is a world-wide trend towards increasing the share of renewable energy 
sources for electricity generation. Since the yield of renewable energy sources like 
e.g. solar and wind power is time-varying and depends on geographic location 
there is an increasing demand for large-scale power transfer over long distances. 
In recent years there have been more and more high-voltage direct current (HVDC) 
links for power transmission commissioned and more are yet to be planned and 
installed. In [1] an extensive overview of the advantages of DC technology for power 
transmission is given. For instance, in Germany there are several embedded HVDC 
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links (i.e. they connect synchronous AC nodes) planned in order to transfer wind 
power from the north to the high load regions in the south, see [2]. 

Optimal power flow is a widely discussed optimization problem to determine 
the optimal operation of controllable equipment in a power grid, see e.g. [3, 4]. 
For transmission grid operation the so-called N — 1 criterion as defined in [5] 
is a requirement. It means that for any single outage of a grid component (e.g. 
a generator, a line or a power transformer) the remaining grid must be able to 
operate safely within its limits. This leads to security-constrained optimal power 
flow problems (SC-OPF). For this class of problems an extensive overview is 
provided in [6]. SC-OPF problems can be formulated either in the preventive way, 
i.e. there is a single set of control variables that has to be feasible in all considered 
contingency cases or in the corrective way, i.e. for each contingency case there may 
be a different set of control variables, and their values can be reached by short-term 
control actions when a contingency case is taking place, see [7]. 

Optimal power flow formulations for hybrid AC/DC grids including security 
constraints are presented e.g. in [1, 8] and [9]. In [10] a corrective SC-OPF 
methodology is employed to investigate the capabilities of embedded HVDC 
systems to compensate for contingency cases in the AC system and thus to reduce 
redispatch by generators. 

In this work we use SC-OPF to study the effects of the N — 1 criterion with 
emphasis on contingency cases within the embedded HVDC systems. To our 
knowledge the effects of N — 1 cases inside the DC systems have not yet been 
studied extensively in the literature. As in [10] we focus on the capabilities of the 
grid components and in particular the embedded HVDC systems to compensate for 
contingency cases and to avoid generator redispatch. 

This work is structured as follows: in Sect. 2 we present the mathematical model 
for SC-OPF for hybrid AC/DC grids. Then, in Sect.3 we describe the case study 
grid and scenario used as well as the implementation of the approach (Sect. 4). The 
case study results are presented and discussed in Sect.5, and we give conclusions 
and suggestions for future work in the final Sect. 6. 


2 The Mathematical Model 


The mathematical model of the SC-OPF is used to perform a steady-state optimiza- 
tion for hybrid AC/DC grids, where the two types of grids are coupled via a set of 
power converters. The model is formulated such that just one single optimization 
problem can be used to solve the AC and DC grids simultaneously; this means in 
particular that no iterative algorithm is needed to reach the final solution. 

The model is written to address one-phase electric systems, i.e. it can evaluate 
a mono-phase system or a multiphase system but in the latter case the system 
must be balanced and symmetric to consider only the positive sequence (Fortescue 
transformation). 
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All the variables and equations are expressed in complex polar form and they 
are differentiable and continuous. This means that discrete variables (e.g. limited 
number of tap changer positions or discontinuous active power range for the 
generators) are not considered. The set of variables are split into two parts: control 
(or independent) variables u representing power setpoints (active and reactive) for 
controllable devices as well as tap changer and phase changer positions for power 
transformers. State (or dependent) set of variables x consisting of voltage magnitude 
and phase at AC nodes, voltage at DC nodes as well as power flows over AC and 
DC branches. 

Regarding the N — 1 criterium, the model is able to take into account the outage 
of a branch (overhead line, power transformer or cable) or a generic equipment 
inside the system (power generator, power converter). Both the preventive and the 
corrective N — 1 methodology can be used. Furthermore, there is no particular 
restriction about the grid topology. This means that it is possible to solve meshed 
AC and DC (e.g. multiterminal DC) synchronous and asynchronous grids. 

The mathematical model can be formulated as the following nonlinear optimiza- 
tion problem: 


Ne 
minimize po f (uo, xo) + I, Pe f (Uc, Xe) 
U0,...UNe 


c=1 
subject to gceluc,x)=0 c=0...N., (1) 
hee, xe) <0 c=0...Ne, 


Au. < Ue — uo < Auc, 


where N. is the number of contingency cases considered (c = 0 corresponds to the 
base case where no contingencies have happened). The set of equality constraints 
&c(Uc, x.) describes the power balance for each node and equipment, while the set 
of inequality constraints h.(u., x.) describes the physical limits of each component. 

The deviation of the control variable in case of corrective N — 1 method can be 
accordingly limited with the upper Au. and lowerAu, bounds (e.g. active power 
setpoint deviation between the contingency cases and the base case). 

The objective function f may be chosen as the generator dispatch cost in 
the system, system losses, voltage magnitude or angle variability, reactive power 
generated or another power system objective. It is also possible to provide flexible 
weights for the different contingency cases under consideration. In this work the 
objective function is defined as the economic dispatch cost in quadratic form: 


Ng 


c=) [oe + BePe + ve P?| 
gl 


(2) 
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Here N, is the number of generators in the grid, P, is the active power generated by 
the generator g, while ag, Bg, Yg are cost coefficients for the individual generators. 
Ap}, c and Ap; , are the active power deviation (positive and negative) for each 
contingency and they are calculated with respect to the base case and weighed with 
the relative linear cost coefficients ct and c. 

The grid model is composed of the following items: nodes (AC and DC), power 
loads, branches and active equipments (generators, converters and reactive power 
sources like e.g. static VAR compensators (SVC)). Two power balance equations 
are given for the AC nodes and one for the DC nodes. In the subsequent paragraphs 
an overview over the models of the grid components is provided. 

The following equation describes the active power balance for an AC node: 


— P(n,cs On,c) + Pg,c + Ph,ac,c — Pn,d = 0, (3) 


where the subscripts n and c denote the n-th node and the c-th contingency case. 
The term p(vn.c, 9n,c) is the active power flow into the n-th AC node where v and 
0 are the voltage magnitude and phase; the term p,,. is the power supplied by the 
g-th active source; the term pp _ac,c is the power supplied by the h-th converter; the 
term pg is the power demand. In the same way it is possible to describe the reactive 
power balance for the AC nodes: 


= g Unes On,c) + 9g,c + In.ac,c — In,d = 0. (4) 


For the DC nodes, only one power balance equation is needed: 


=> P(Un,c) + Ph,dc,c = 0, (5) 


since in this case no power generators and loads are present. For both AC and DC 
nodes the voltage magnitude is constrained to its lower and upper limits. 

The branches are modeled as generic passive components able to transfer current 
from one node to another according to Kirchhoff’s circuit laws. The double bi-pole 
model with PI scheme has been used. This model can be applied to overhead lines, 
cables, power transformers (with and without tap changer and/or phase shifter) 
and to DC overhead lines and cables. Without loss of generality we assume the 
tap changer and phase shifter to be installed on the from-node of the branch and 
constrained with their limits. Both can be modelled with a pure gain 6 and a complex 
rotation e/®. For a generic branch b, contingency c and defining i as from-node and j 
as to-node, it is possible to describe the relationship between the current and voltage 
phasors as follows: 


is} | Wi Vance! | | vj | _ 6 
ij Yiröp,.ce I< Yy vj 
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Exploiting the symmetry of the model, the matrix component with subscript ii 
is equal to the component with jj and, in the same way, the component with ij 
is equal to the component with ji. Those values arise from the longitudinal and 
transversal susceptance and conductance parameters. The equations that describe 
a double bi-pole model are used to model the active and reactive power flow 
through the branches, and thus those flows are used inside the nodes power balance 
equations (3), (4), (5). Obviously, the power flow (active power or apparent power or 
current module) can be constrained at both sides of the branch with an upper limit. 

The power generators are generic active components able to inject active and/or 
reactive power into the system. The generator capability is defined as a domain 
created by the intersection of a rectangle (active and reactive upper and lower 
bounds) and a circle (apparent power bound). Additional linear bounds between 
active and reactive power can be defined in order to customize the capability domain. 

In a similar way, by neglecting the active power, voltage controlling reactive 
power sources can be used inside the model: in this case, a linear relationship 
between voltage magnitude and reactive power output is used. 

The power converters are able to transfer power in both directions between AC 
and DC nodes. The converter model takes into account the power balance between 
the AC and DC ports with the losses due to its power electronics. Moreover, the 
converters may have capability domains described in the same way as for generators. 
The converter power balance equation is written as follows: 


Ph,ac,c + Ph,dc,c + Ph,loss,c = 0, (7) 


where the subscripts loss denotes the power losses inside the h-th converter and the 
c-th contingency. The AC side current is determined by the following equation that 
describes the power balance in the AC converter side: 


fie, 2 
Pn,ac,c + An,c 


J3 Un,ac,c 


Finally, the power losses are calculated by a formula that creates a quadratic 
relationship between the current iy ac,c and losses pn joss,c (see e.g. [11] and [12]) 
with parameters œ}, Bn and yp that may be individual for each converter: 


= th,ac,c = 0. (8) 


. 2 
an + Bnin,ac,c a Yhln,ac,e — Ph,loss,c = 0. (9) 


This model is suited for steady-state analyses, so that transient behaviour, e.g. 
fault and stability analysis cannot be investigated with it. 
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3 Case Study 


The case study in this work is based on the IEEE1 18 test grid (see [13]). 

The generator costs have been modified in order to encourage power transfer 
from the (low cost) western zone to the (high cost) eastern zone (see Fig. 1). In 
particular, the generation cost in the western zone have been set to zero to emulate 
renewable generation, while the cost coefficients in the western zone have not been 
modified. 

The focus of our analysis lies on one additional power corridor between the nodes 
8 and 65 (highlighted in red in Fig. 1). There are three options for this considered 
corridor: an AC corridor, a monopolar HVDC corridor as well as a bipolar HVDC 
corridor with identical total power ratings of 500 MVA and total length of 500 km. 
The electrical parameters of the AC corridor as well as the resistance of the DC lines 
are taken from [14]. Due to the huge length of the corridor the AC line has been split 
in three line segments. The HVDC converters exhibit a loss model with a quadratic 
relationship to the AC current (see [15]). 

For all three configurations the effects of security constraints on the corridor 
are analysed. For both the AC corridor and the monopolar HVDC corridor the 
considered contingency implies the outage of the whole corridor, while for the 
bipolar configuration an N — 1 contingency case involved the outage of one of the 
two poles. 

In order to discourage redispatch of generators in the contingency case, we 
require the generator active and reactive power output to be identical in the base 
case and the contingency case, while HVDC converter setpoints may differ between 
base case and contingency case. This means that the generators follow a preventive 
security strategy while HVDC converters follow a corrective security strategy. The 
same behaviour can be achieved by imposing very high redispatch cost on the 
generators. The possibly changing losses after an outage are balanced by adapted 
voltage profiles. In addition, the setpoints for transformers are fixed. 


4 Implementation 


A software tool for solving the security-constrained optimal power flow problem as 
described in Sect. 2 has been developed in C++ using Microsoft® Visual Studio™. 
Inside this software the open-source library IPOPT (see [16]) is used as solver 
for non-linear programming based on an interior point method. To speed up the 
optimization process and to increase accuracy, we derived the analytic forms of 
Jacobian and Hessian matrices of the objective function and the constraints and 
pass them directly to the solver instead of using numerical differentiation. 

The input data as well as the optimization results including all control and 
state variables are given as text files. In the results discussion below (Sect.5) we 
summarize these results in the form of the key performance indicators dispatch 
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cost, power flow over the considered corridor, active power losses and voltage angle 
deviation. 


5 Results 


The result of the SC-OPF is analysed in terms of dispatch cost (Fig. 2), power flow 
over the considered corridor (Fig. 3), active power losses in the system (Fig. 4) and 
voltage angle difference along the corridor (Fig. 5). 

The first observation is that in the cases without the N — 1 criterion, all the values 
for the monopolar and the bipolar HVDC corridors coincide, which is obvious since 
they have identical parameters and rating. When looking at the dispatch cost in 
Fig. 2 it can be observed that the cost is lowest for both the HVDC variants without 
N — 1 criterion. With N — 1 criterion both HVDC variants have lower dispatch cost 
than the AC variant, while the bipolar configuration has the lowest dispatch cost. 
This can be explained by the higher efficiency of the HVDC systems compared to 
the AC corridor (lower losses over long distances) and the higher redundancy of the 
bipolar configuration compared to the monopolar one. This way more power can 


Economic Dispatch 
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Fig. 2 Dispatch cost relative to bipolar HVDC scenario without N — 1. Red, dark blue: AC line 
with/without N — 1; violet, yellow: monopolar HVDC with/without N — 1; light blue, green: bipolar 
HVDC with/without N — 1 
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Fig. 3 Power flow over corridor relative to bipolar HVDC scenario without N — 1 (457 MW) 


be transferred from the low cost generation zone to the high cost generation zone 
which enables a generation dispatch at lower total cost. 

This explanation is confirmed by the power flow over the corridor (Fig. 3), where 
the bipolar HVDC configuration is able to transfer the most power from the low cost 
generation zone to the high cost generation zone. The monopolar HVDC corridor 
is transferring significantly less power in the N — 1 case since it has to prepare 
for the outage of the full corridor. The AC corridor is not even able to maintain the 
power flow direction in the N — 1 case. The significant deviations in the power flows 
between the AC and monopolar HVDC configuration can be explained by voltage 
limits and differences in the generator dispatch. 

The active power losses in the system can be seen in Fig. 4. There we did not 
observe a clear pattern, which can be explained by the fact that in all the scenarios 
the generator dispatch and thus the power flows through the whole system are very 
different, and the losses were not the objective of the optimization performed. 

When looking at the angle deviations in Fig.5 is can be observed that the AC 
corridor exhibits a large angle deviation in the case without security criterion, while 
all the HVDC configurations are able to maintain smaller angle deviations which 
indicates a higher system stability. 

In Table 1 the voltage magnitudes and reactive power injected at the corridor 
terminal nodes 8 and 65 are listed for the N — 1 scenario before and after the 
contingencies. It can be observed that the bipolar HVDC configuration is able to 
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Fig. 4 Active power losses in the system relative to bipolar HVDC scenario without N — 1 
(120 MW) 
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Fig. 5 Angle Deviation over corridor in degrees 
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ae 1 Voltage magnitude ve) | V(65) | Q(8) Q(65) 
at modes 8 al €5 forallihre Case [pu] [ipu] [MVA] | vag 
corridor configurations in the AC; base 0.980 | 0.968 98.2 118.6 
base case and the contingency AC, cont. 1.046 | 0.965 0.0 0.0 
case of the N — 1 scenarios DC MP; base | 1.039 | 1.012 |-19.0 | —96.1 
DC MP, cont. | 0.973 | 1.021 0.0 0.0 
DC BP; base | 1.004 | 1.023 | —58.7 —5.1 
DC BP, cont. | 0.998 | 1.023 | —56.5 0.0 


maintain a stable voltage level close to the nominal voltage by providing reactive 
power to the grid. In contrast, for the AC configuration the voltage magnitude 
at node 8 changes significantly because of high generation in the region. The 
monopolar HVDC configuration is also stabilizing the voltage by reactive power 
injection but in the outage case there is the same effect as for the AC corridor. 


6 Conclusion and Future Work 


A mathematical model and a software for security-constrained optimal power flow 
for generic hybrid AC/DC grids has been developed. In this work it has been used 
to study the effects of the N — 1 criterion on embedded HVDC systems. 

The above results show that for the bipolar HVDC configuration the power plant 
dispatch with lowest cost could be established when the N — 1 criterion is employed. 
This study has been performed on the IEEE118 grid with a single corridor under 
consideration. The results indicate the benefits of employing bipolar HVDC systems 
in order to increase the social welfare. 

In order to draw conclusions for real transmission systems, additional studies on 
real grid configurations with possibly several embedded HVDC links would have to 
be performed. 

Additionally, an investigation including contingency cases both in the AC and 
DC parts of the power system would be interesting to investigate the capabilities 
of embedded HVDC to compensate for general N — 1 cases and thus to reduce 
redispatch cost for grid operators. 
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1 Introduction 


The increasing penetration of Renewable Energy Sources (RES) leads to a power 
system operation closer to its operational limits. Enhanced methods and tools will be 
crucial for a secure, but also efficient and cheap planning and operation of the power 
grid. A Transmission System Operator (TSO) will have larger assets of controllable 
devices at hand, such as Voltage Source Converter (VSC-) based HVDC-systems 
and energy storage systems. Combined with the need to assure N-1 security, this 
requires a Security-Constrained Dynamic Optimal Power Flow (SC-D-OPF), which 
incorporates both N-1 security and multiple time steps into the optimization. First 
research in that area was done in [1], where not only power but also necessary 
energy reserves can be determined to ensure N-1 security. In [2], successive linear 
algorithms and approximated power flows are used to solve large-scale SC-D-OPF 
problems in a European context. The authors of [3] propose an SC-D-OPF model 
including uncertainty of wind power or equipment availability, which is solved in 
a two-stage stochastic program. In [4], SC-D-OPF is used as inner iteration in a 
hierarchical approach to optimize a central deployment signal which is sent to the 
units able to provide a reserve. Furthermore, [5-7] include energy storage systems 
in an SC-D-OPF calculation. A further extension is shown in [8], where SC-D-OPF 
is solved in a hybrid AC-DC grid with energy storage. 

Due to the increasing complexity of the power system and an operation closer to 
network limitations, central coordination in large scale networks comes with major 
computational burdens. Furthermore, privacy can be a concern for each transmission 
system operator (TSO) controlling a certain region. Subsequently, the interest in 
distributed optimization, also referred to as multi-area optimization has substantially 
grown in recent years. An early overview of distributed OPF algorithms can be 
found in [9] and the most recent developments are examined in detail in [10]. The 
Alternating Direction of Multipliers Method (ADMM) [11], has become a popular 
branch to tackle the non-convex AC OPF problem [12-14]. Each area uses variable 
duplicates from neighboring areas and is solved to optimality. Penalty terms, which 
are calculated from the coupling variable deviation and Lagrangian multipliers, 
are added to the objective function to push two neighbors towards a common 
boundary solution. ADMM was applied to a hybrid AC-DC grid in [15]. Multi- 
area optimization is applied to D-OPF with storage in [16] using OCD, but to the 
best of our knowledge, no attempts have been made toward distributing SC-D-OPF. 


2 Security-Constrained Dynamic Optimal Power Flow 
(SC-D-OPF) 


The optimization horizon is defined by a set of time steps 7 = {1,..., T} with 
T the total number of considered time steps. Furthermore, we define a scenario set 
C = {0,..., C}. Here, scenario c = 0 defines the base case, i.e. the original fault- 
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free network. A total of C possible contingencies is added by including modified 
versions of the base case to the scenario set. A modification could be the outage of 
an AC line, a DC line or a converter. With x" the optimization variables of time 
step ¢ and scenario c, the full set of optimization variables is collected in 


SG ET ENT a): (1) 


The full problem has the form (2): 


eee tic tc 

minimize 22 p fS) (2a) 

subject to hpase(x"®) <0 YteT,cecC (2b) 
hiya” Heil (2c) 
hyi (x, x") <0 WteT,cEC\0. (2d) 


It minimizes the sum of weighted costs f over all scenarios. Factor p! defines 
a probability of each contingency to enable a risk-based OPF. Constraints /pase 
describe the basic constraints of a conventional OPF, which must be fulfilled during 
each scenario (t, c). That includes AC and DC node power balance; thermal AC 
and DC branch limits; AC and DC voltage limits; quadratic loss function of AC- 
DC converters; and operational limits of generators, converters or sheddable loads. 
Time-dependent constraints, such as storage energy limits or generator ramping 
limits, are collected in hayn. Those constraints only apply to the base case (c = 0), 
since set points after an outage are coupled solely to the respective base case of 
the identical time step. This N-1 coupling is modeled by the constraints hn-ı and 
guarantees N-1 security. Here, the coupling between each contingency and the base 
case is explicitly modeled. That is, set point deviation limits from before to after 
the occurrence of an outage can be defined. Note that for the sake of readability, we 
omit equality constraints which can be interpreted as reformulated inequalities as 
well. We refer to [8] for more details on the modeling. 


3 Distributed Formulation 


For the sake of readability, we write (2) in the form (3): 
minimize f(x) (3a) 
x 


subject to h(x) <0. (3b) 
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Let R non-overlapping regions be defined in R = {1,..., R}, that is, each node 
belongs to exactly one region. An equivalent problem as (3), but in separable form, 
can be written as (4): 


Er j 

minimize XO fa) (4a) 
keR 

subject to har) <0, WkER (4b) 
X A=. (4c) 
keR 


Here, x € R*, hg : Rk > R™ and Pe: Rk — R are optimization 
variables, non-linear constraints and objective function, respectively, in a region 
k. Matrix Ak € R"*’« maps x, onto the full set of n coupling constraints and 
enforces consensus between regions. Thus, (4b) forms independent local SC-D- 
OPF constraints and (4c) form the coupling constraints which recuperate a feasible 
overall power flow. 


4 Network Decomposition in AC-DC Grids 


To allow for a separable formulation, the network is first partitioned and then 
decomposed. 

Network partitioning can be crucial for distributed algorithms to achieve good 
performance. We believe that network partitions are inherently given by structural 
responsibilities. For example, a region or control area could represent one TSO, 
multiple TSOs or a whole country. Thus, the number and dimension of AC regions 
are historically known. However, responsibilities in overlaying DC networks are 
yet to be defined and various options are thinkable [15]. We choose a Shared-DC 
approach. Here, we define R = RA© hybrid AC-DC regions, where each AC region 
may also contain DC nodes. Thus, each existing TSO gains control over converters 
in its own control area. Let Mg identify all nodes in region k. Herewith, N collects 
all AC nodes and NPE all DC nodes in region k. 


4.1 Decoupling of AC or DC Tie Line 


In principle, the decoupling of AC and DC tie line is identical, except for extended 
consensus constraints due to complex voltage and power in the AC case (Fig. 1). 
The original line is cut into two halves; auxiliary nodes (m,n) and auxiliary 
generators (a, b) are added at both open ends. Thus, node sets are augmented to 
Na = {Na,m}, Ng + {NpB, n}. To guarantee a feasible power flow, the voltage 
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Fig. 1 Decoupling model of a tie line (AC or DC) between nodes i and j. Two auxiliary nodes (m 
and n) and two auxiliary generators (a and b) are added at the middle. Reactive power source Qg 
is only added for an AC tie line 


must be equal at nodes m and n. Furthermore, the generators must produce the same 
amount of power of the opposite sign. In the case of an AC tie line (i € M re Je 
N ok this leads to boundary conditions including complex voltage and both active 
and reactive power. For time step and scenario (f, c), we have: 


Vato! = [Vaca (5a) 
LViem = LV ncn (5b) 
Poa = Po (Se) 
Oga =: (Sd) 


In the case of a DC tie line (i € NRS, TEN, - only real voltage and active power 
must meet the constraints. For time step and scenario (f, c), we have: 


të t,c 
Voom = Yocn (6a) 
DE _. HE 
Pga = Por (6b) 


Equations (5) and (6) are the only boundary constraints. 


4.2 Building Consensus Constraint Matrix A 


Since oe uses the augmented node sets of region k, the auxiliary generator variables 
are included. Thus, consensus constraints (5)-(6) between region A and B can be 
written to 


Alex’ + A xg = 0. (7) 
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Fig. 2 Schematic problem decomposition for two areas (blue and green) under consideration of 
three time steps and two contingencies. A dashed box represents an optimization problem. (a) 


Centralized SC-D-OPF. (b) Distributed SC-D-OPF 


This must be fulfilled for every considered time or contingency scenario (Fig. 2). 


Then we have 


In a more compact form, (8) can be written as 


= Arxk=0. 


ke{A,B} 


5 Implemented ADMM Algorithm 


(8) 


(9) 


The general idea of ADMM, see also [11], is the following. Augmented regional 
OPFs are solved and the deviation of boundary variables from a fixed auxiliary 
variable z, which is information stemming from neighboring regions, is penalized. 
The regions then exchange information and z is updated. The update is re-distributed 
to the local agents (areas) for a new OPF calculation until consensus between 
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regions is achieved. It is constructed an augmented Lagrangian of the form 


LS Inn + Ag Akk — Zk) 
keR 


+ Bias — ally} (10) 


with penalty parameter 0 € R. Dual variables of the consensus constraints are 
denoted with A, € R"*!, and W € R”*” is a positive definite, diagonal weighting 
matrix,! where each entry is related to one coupling constraint. We refer to W(S) 
if the entry is related to a power variable, and to W(V) if the entry is related to 
a voltage variable. In ADMM literature, W is the identity matrix. The main steps 
during one iteration of the solving process are 


l. x = argmin ex L(x, z, A) (11a) 

2. z = argmin ez L(x, Zz, A) (11b) 

3. à <—A+ PWA(x — Zz) (llc) 
with z = [z] set Zh) The first step (11a) minimizes a non-linear problem, 
where constraint region 

x= {x |g (xx) <0,Vk € R} (12) 


enforces local constraints (4b). Since z is fixed, (lla) is in fact a series of 
R independent problems which can be calculated in parallel. The second step 
minimizes a coupled quadratic problem, where 


Z= 42| >) Aree =0 (13) 
keR 
enforces consensus of auxiliary variables z. In fact, (11b) calculates the average 


value between two consensus variables of neighboring regions [11, 13] and the 
problem can be reduced to 


p 
z = argminez ) | Z|lArt — zoll: (14) 
keR 


IA weighted norm is calculated with ||X | 2, = X'Wx. 
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z sp 
G2 

Qa 

Fig. 3 Test system with embedded HVDC, energy storage and RES. Partitioning into 3 control 

areas 


The weighting factors o and W can be neglected in (14), if the same weight is 
assigned to two coupled variables, which is a reasonable choice. Once a region has 
gathered necessary neighbor information, this step can be calculated locally as well 
[12]. In the third step (1 1c), dual variables are updated based on a weighted distance 
between x and z. Again, with given local xx and zx, each region can calculate Ax 
independently (Fig. 3). 

In ADMM, an update rule on p can be useful to enforce consensus. A pp € R 
is assigned to each region, which can be increased depending on the local residual, 
see (18b). If the residual has not decreased sufficiently compared to the previous 
iteration (indicator 0 < © € R < 1), the penalty is increased by a constant factor of 
t € R > 1. The penalty parameter must be chosen carefully since it is widely known 
to be crucial for good convergence behavior [17]. An overview of the implemented 
ADMM is given in Algorithm 1. 


6 Results 


6.1 Test Scenario 


We use the illustrative AC-DC test system from [15], see Tables | and 2 for line and 
generator parameters, respectively. The generator at Node 5 is interpreted as a wind 
park and a PV park is connected to Node 3. RES upper power limit and load are 
variable over time. Forecast power profiles are taken from the Belgian transmission 
system operator Elia, which provides detailed RES generation and load data on its 
website. Adapted profiles for usage in the 5-bus test system are shown in Figs. 6 and 
7. We assume quadratic generator cost functions as in [18] (see Table 3) and one 
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Algorithm 1 ADMM 


1: Initialization: Weighting matrix W, tolerance e€; for all k € R: initial guesses zg, penalty 
parameters ok = p, dual variables àg = 0, local solutions x; = 00, local residues T% = 00. 

2: while ||Ax||.. > € and ||x — z||o > € do 

3: Solve for all k € R the decoupled NLPs 


minimize fe (e) + Ag Auıı + PE IIA 01%, (15a) 
subject to hg (xx) < 0 (15b) 


4: Solve the coupled averaging step 


minimize > ||Ax (xe — zoll (16a) 
É keR 

subject to D Arz =0 (16b) 
keR 


5: Update dual variables for all k € R 


Ak — Ak + peWARK(XK — zk) (17) 
6: Calculate local residues and penalty parameter updates for all k € R 
TE = lA — ze) I loo (18a) 


if T} < en 
pe ee eee (186) 
Tpk otherwise 


7: Update I, = I: forall k ER 
8: end while 


cost coefficient for all reactive power injections. RES power injection is prioritized 
by assigning a low fuel price. Two energy storage systems with characteristics from 
Table 4 are connected to Node 3 and 5, respectively. Line flow limits of 400 MVA 
for Line 1-2 and 240 MVA for Line 4-5, respectively, are used. Furthermore, the 
capacity of each VSC is set to 200 MVA. Allowed voltage ranges lie between 0.9 
and 1.1 pu on the AC side, and between 0.95 and 1.05 pu on the DC side. 


6.2 Multi-area SC-D-OPF 


To show applicability of the distributed approach to N-1 secure and dynamic 
OPF problems, a small example is presented in the following. An SC-D-OPF 
dispatch optimization is performed for the scenario set C = {0, Line 1-2} and 
T = {1,..., 4}. Storages and VSCs are allowed curative control, that is, operating 
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Fig. 4 Convergence behavior of distributed SC-D-OPF in the 5-bus system with 4 time steps and 
1 contingency. Left: consensus error, right: cost error. Both relative to central solution 


set points are allowed to move freely to a new set point after an outage. Contrarily, 
conventional generators are operated preventively, that is, the power set point must 
be valid both before and after the outage. However, to cope with changing system 
losses after an outage, a deviation of | % of installed capacity is allowed. 

We use ADMM settings p = 500, t = 1.02, © = 0.99, W(S) = 1 and W(V) = 
100. 

General convergence behavior is shown in Fig. 4. Consensus error and objective 
suboptimality show similar behavior compared to conventional OPF cases. Consen- 
sus error (|| Ax ||o0) which depicts cross-border feasibility, falls below the threshold 
of 1074 after 344 iterations. The deviation from centrally computed costs f* is 
denoted with f = |1 — f/f*|. The objective value error falls below 0.01 % when 
convergence is achieved. 

The progress of storage variables is shown in Fig.5 for all time steps after the 
outage. Dashed lines denote results from the centralized solution. Storage systems 
provide positive (ESS 2) and negative (ESS 1) power reserve after the outage. The 
amount of reserve is increased with time for a growing network stress level, but 
energy levels remain within limits. The distributed solution is represented by solid 
lines. Those oscillate around and eventually approach the target levels. 


7 Conclusion 


This paper presents an approach to calculate SC-D-OPF in a distributed manner. 
That is, each control area computes a local N-1 secure and dynamic optimal dispatch 
which respects boundary constraints from neighboring areas for each time step and 
topology. After iteratively exchanging information and updating the local solution, 
a consensus is found which allows for a feasible cross-border power flow for all 
possible scenarios. Additionally, the result is close to the central optimal solution. 
Optimal curative storage or VSC set points can be determined for a contingency 
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Fig. 5 Convergence of storage power and energy levels after the outage 
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scenario which is in a foreign area. In this work, we apply the method to a small 
test system — the next steps are to increase system size, time horizon and number of 
contingencies. 


Appendix 


Table 1 Line parameters Line! R fpu] |X [pu 
with base power 100 MVA 


0.00304 | 0.0304 | 0.00658 


0.00064 | 0.0064 
0.00108 | 0.0108 


0.00297 | 0.0297 | 0.00674 


0.00297 | 0.0297 


Table 2 Generator Goce oP 
parameters. (Pg, Pg) in lOc | 
[MW] and (Qa, Qa) in a fo [170 | -127.5 | 
[Mva] 3 |o [520 | -390 


4 [o 1200 |-150 


Table 3 Generator cost Cost coefficients 


coefficients Eea 
Generator ba [ uw £] 
G1 0.010 15 
G2 0.011 30 


G3 40 


Table 4 Parameters for each 
energy storage system 


200 MWh 
100 MW 
20 MWh 
180 MWh 
95% 


Capacity 
Maximal power 
Minimal energy 
Maximal energy 
Efficiency 
Costs 

— discharge 10 €/MWh 


— charge 0 €/MWh 
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Fig. 6 Load profiles 
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Fig. 7 RES profiles 
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Problems in Parallel 
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Thomas Leibfried, Wolf Fichtner, Valentin Bertsch, and Vincent Heuveline 


Abstract We propose a parallel solver for linear systems of equations arising 
from the application of Primal Dual Interior Point methods to Dynamic Optimal 
Power Flow problems. Our solver is based on the Generalised Minimal Residual 
method in combination with an additive Schwarz domain decomposition method 
as preconditioner. This preconditioner exploits the structure of Dynamic Optimal 
Power Flow problems which, after linearization, is given as a matrix with large 
diagonal blocks and only a few off-diagonal elements. These elements correspond 
to intertemporal couplings due to ramping and energy storage constraints and are 
partially neglected in order to solve the problem in parallel. We test our method on a 
large-scale optimisation problem and show that a parallel speedup can be obtained. 
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1 Introduction 


While today’s power grid infrastructure has been designed for centralised power 
generation in conventional power plants, the ever increasing share of renewable 
energy sources (RES) leads to an increasingly uncertain, volatile and decentralised 
generation. In order to ensure a reliable grid operation and to maintain today’s 
security of supply, it is necessary to develop methods to simulate accurately a power 
grid at high temporal resolutions. This can be done by efficient methods for power 
grid optimisation, including an accurate consideration of the non-linear and non- 
convex AC power flow constraints. And these are to be improved. 

The task to find the optimal operating state of a given power grid, also known as 
Optimal Power Flow (OPF), is formalized as the minimisation of a cost function 
with respect to a vector of continuous optimisation variables like the generated 
active power and the generated reactive power. Dynamic Optimal Power Flow 
(DOPF) considers several OPF problems, each corresponding to one specific 
point in time. In this case, the power grid’s power demand is time-dependent and 
additional constraints must be taken into account. These constraints couple some 
of the optimisation variables corresponding to different time steps. Among others, 
these couplings are introduced by energy storage facilities and ramping constraints 
for conventional power plants [20]. Since DOPF is a large-scale non-linear and non- 
convex optimisation problem, solving it in parallel is of crucial importance. 

For this kind of problems Primal Dual Interior Point Methods (PDIPM) have 
proven to be among the most powerful algorithms. The main computational effort, 
when PDIPM is applied, lies in solving linear systems of equations [14]. And this a 
task suitable to be carried out in parallel on multi-core CPUs. 

There exist a few works on solving linear systems, that arise from PDIPM 
applied to DOPF, in parallel. One approach is based on parallel matrix factorization 
by means of Schur complement techniques [6, 20]. Other strategies use Benders 
Decomposition to decompose the complete optimisation problem into smaller ones 
[1]. 

Our contribution is the use of a Schwarz Domain Decomposition Method as 
preconditioner for Krylov-type iterative methods for solving these linear systems 
of equations in parallel. The original Schwarz method was formulated in 1870 as 
theoretical tool for proving existence of elliptic partial differential equations (PDEs) 
on complicated domains [17]. Later on, modifications of it have been used as stand- 
alone iterative methods for solving PDEs and have become a standard technique for 
preconditioning Krylov-type methods in the context of PDEs [18]. We apply this 
technique to DOPF by decomposing the time domain of interest into several smaller 
subdomains, which allows the use of multiple processors for computation. 
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2 Dynamic Optimal Power Flow 


We will now formulate the DOPF problem in a mathematically rigorous way. This 
formulation has already been stated in more detail by several authors, e.g. in [12, 20] 
and [16]. We would like to mention that we consider a continuous formulation of 
the DOPF problem, i.e. one without any discrete decision variables. A non-linear 
mixed-integer formulation of an OPF is considered, for example, in [19]. 


2.1 Formulation of the DOPF 


In order to formulate the DOPF problem for a given time period of interest [0, 7], 
we consider a uniform partition 0 = Ti < T2 <... < Ty, = T with constant step 
size T = T; — T;_; for allt € T \ {1}, where T := {1,..., Nr}. 

The power grid consisting of Ng nodes denoted by B := {1,..., Ng} and Ng 
transmission lines denoted by E C B x B, is described by an admittance matrix 


Y=G+jBeCs*%s 
(1) 
with G, Be RNB*4B and j = /—1. 


The admittance matrix is symmetric, i.e. Y = Y T and furthermore Yır = Yu + Oif 
and only if there is a branch connecting node k and /, i.e., (k,/) € E and (l, k) € E. 
Therefore, Y is a sparse matrix for most real world power grids. 

The complex voltage at node k € B for time step r € T is given as the complex 
number 


Vi =E +jF! (2) 
with real part Ef € R and imaginary part F! € R. We define E’ := [E}]keg and 


F' := [Fi)keB- 
Moreover, we define the following sets, which number: 


C := {1,..., Nc} conventional power plants 
R := {1,..., Nr} renewable energy sources 
S := {1,..., Ns} storage facilities 


LS := {1,..., Nis} load shedding possibilities 
DC := {1,..., Npc} generators modelling DC lines 
V := {1,..., Ny} generators modelling import/export 
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The resulting vector of variables corresponding to active power and energy for 
time step f is given by 


[Pe ili eC (injected conventional power) 
[P ri Jier (injected renewable power) 
i g,ilies (injected storage power) 
[Pi ilies (extracted storage power) 
pi _ | [SOCIlics | (stored energy) N 
[Pi liets (load sheds) 
[Piep ilieDC (injected power by DC lines) 
[ Pic il iepc | (extracted power by DC lines) 
[ Ph ev (power export) 
[P!, „liev (power import) 
with P’ ER". 
In a similar way, the vector of reactive power injections is defined by 
[0% ‚liec 
(a! ‘lier 
Q= | 105, ‚lies | € RN. (4) 
LO, ilies 
[Q} ilies 


Every variable related to active and reactive power can be assigned to a specific node 
of the power grid by means of the power injection and extraction matrices Cp € 
{—1, 0, 1}%8*? and Cg e {—1, 0, 1}%8*%e [21, p.23]. Since the states of charge 
(SOC) of storage facilities do not directly influence the power flow equations, the 
corresponding rows in Cp are equal to the zero vector. 


2.1.1 Constraints 


At first, we consider the equality constraints. We begin with the power flow 
equations: Denoting the active and reactive power load for time step t by Li and 
Lie respectively, and assigning the loads to the nodes with the help of the matrix 


Cz € {0, 1}Y8*NL, the AC power flow equations are given by [22] 


CpP' - CLL) — Pfp(E', F') = 0 e RS, (5) 


CoQ‘ — CLL) — Pf,(E', F') = 0 ER, (6) 
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where 


Ng 
Pfp (E, F) = > Ex (Gu Er — Bu Fi) + Fr(Grı Fi + Bu Ei), 
[=1 


NB 
Pfq,x(E, F) =) Fx(Gu Et — Bu Fi) — Ex(Gu Fi + Bu Er), 
l=1 


which describe the sum of power flows at node k. Note that we dropped the index t 
to improve readability. 

In the following, we abbreviate the equations (5) and (6) with the help of the 
following “symbolic” equation: 


AC(E", F', P', Q') =0. (7) 


Besides the power flow equations (7), one has to take into account an equality 
constraint for the slack bus, given by F i = 0 for all ¢, and for generators k,l € DC 


modelling a DC line Kl: 
Den = —vpc Pie | and Pier = —vpc Pick : (8) 


Here, line losses are taken into account by means of the factor vpc € (0, 1). 
Moreover, an additional set of equality constraints has to be imposed for 
modelling the state of charge, SOC!, of storage facilities. These constraints are 
given by 
SOC} = SOC! + tv P! 


= 
apo, E 


eae forallieS (9) 


with factors vj, vg € (0, 1) describing the efficiency of the loading process P! qi and 
generation process P! gi respectively. 

Secondly, we take a look at the inequality constraints. Due to technical restric- 
tions arising from power plants, renewable energy sources, storage facilities and DC 
transmission lines, it’s necessary to impose lower and upper bounds for active and 
reactive power injections: 


er = P' < Plax and Qmin 2 Q' < Qmax . (10) 


In our application the only bounds, which are time dependent, correspond to the 
power injected by renewable energy sources, i.e. Pt. Also the node voltages and the 
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active power flow over transmission lines, denoted by pı, must be constrained: 


i (Py ee ye, (11) 
pei(Ej, Ej, Fi, Ff) < Pfmax,x for all (k, 1) € E. (12) 


Finally, we impose ramping constraints for conventional power plants to bound the 
rate of change for injected power between consecutive time steps by 


|P; Pil <tP!_,a forallt eT \ {1}. (13) 


a is the ramping factor in percentage per unit time. 


2.1.2 Cost Function 


The cost function f accounts for expenditure and income related to active power 
processes for the complete time interval [0, 7] and is composed of contributions 
from each time step t € T [5]: 


Nr 


FP) =Y fP) (14) 


t=1 


with 


AP) =) _(ac,.2(Pi )” + acia Pt) 


ieC 

+ X (aria (PL? + aRiiP';) 
icR 

+F X (aex, ia PhD =F aex i1 Pixi) (15) 
ieV 


2 
+ Y nal) + tim, i1 Pimi) 
ieV 


F > (ats,i,2(P};)° + axs,i,1 Pfi) 
iELS 


and coefficients ax, x, € R. 
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2.1.3 Optimisation Problem 


With the definitions of the previous sections at hand, the DOPF problem can be 
stated in an abstract way: 


min F(x) st. 
eG) = 0, for allt € T 
gga, x 1) =0, for allt € T \ {1} (16) 
ha) < 0, for allt € T 
Woe a YS 0, for allt € T \ {1} 


where the vector of optimisation variables is given by 
x=([x'her € R” and x' = (£ Fi pt o') ER, (17) 


Moreover, we used and will use the following definitions: n*’ := 2Ng + Np + N Q 
and n* := Nrn* and, as of now, let nò’ be defined as the number of equality 
constraints corresponding to time step t, i.e. A’ € R^" with A’ = (At, a) denoting 
the Lagrangian multipliers corresponding to gi and Bes respectively. Analogously, 
u = (uh, uR) € R”"" denotes the Lagrangian multipliers for the inequality 
constraints hi, and hp. Consequently, we define n* := Nrn*" andn” := Npn”. 

Note that the ramping constraints hr given by (13) and the storage constraints 
gs given by (9) establish the only couplings between the variables corresponding to 
different time steps. Without them, a solution to (16) could be obtained by solving 
Nr independent OPF problems. The optimisation problem (16) is non-linear due 
to the AC equations and due to the voltage and line flow inequality constraints.! 
Moreover, the latter ones and the AC equations are the sources of non-convexity 
in (16). 

To simplify the notation in the following sections, we abbreviate the optimisation 
problem (16) to 


min f(x) st g(x) =0, h(x) <0 (18) 
x 
with quadratic functions 


f:R" > R, g: RË > R”, h: R” > R". (19) 


'We have not explicitly stated an equation for the power flow pxı through the transmission line 
(k, I). Please refer to [13]. 
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3 Primal Dual Interior Point Method for DOPF 


In this section, we briefly describe the application of the Primal Dual Interior 
Point Method (PDIPM) to a DOPF problem. The description is in line with the 
implementation found in the MATLAB package MATPOWER [21]. 


3.1 The PDIPM Algorithm 


For the general optimisation problem (18), the corresponding Lagrangian function 
is 


L: R? xR” xR” SR m 
(x, à, U) > f(x) +AT g(x) +u ha). 


Assuming additional constraint qualifications, e.g., the linear independence 
constraint qualification (LICQ) [14], for every local minimum x* of (18), there exist 
corresponding Lagrangian multipliers A*, u* such that (x*, A*, *) are in accord 
with the Karush-Kuhn-Tucker (KKT) conditions [5]: 


Vx L(x, A, u) 
g(x) 
Fo(x, s, à, u) := =0, s,u 2 0, 21 
o(x, s, À, u) Avs S, u (21) 
Su 
where s € R”“ is a vector containing slack variables and S € R""*”" denotes a 


diagonal matrix whose diagonal elements are Skk = sx. 

The PDIPM is an algorithm to find solutions to Eq. (21). Mathematically it solves 
a slightly modified version of this equation by applying to it Newton’s method. 
The modification consists in adding a “term”, which decreases with iteration index 
k. This term keeps the slack variables s away from zero and, thus, the inequality 
constraints inactive. This means that solutions (x ©, s®, m, u®) to this modified 
version of Eq. (21) are inside the feasible region and not on its surface, hence, the 
method is called interior point method. More rigorously formulated, we find 


T 
E, (x, s, à, 1) = Fo, s, à, W) — (000 ye) =0, a 


s,u = 0, 


A Domain Decomposition Approach to Solve Dynamic Optimal Power Flow... 49 


for a sequence of barrier parameters 


(uR-DyT g(k-D 
Yk := 0 Fin 
where o = 0.1 is the centering parameter and (u“—))? s*— is called comple- 
mentary gap. Additionally, e = (1,..., 1) € R"". For details we refer the inclined 
reader to [5, 8] and [21]. 


Remark A very basic PDIPM is applied, i.e. there is no globalisation strategy. 
Since we focus on parallel linear algebra techniques for DOPF problems, a locally 
convergent algorithm is sufficient to our purpose. 


3.2 KKT Matrix in PDIPM 


The application of Newton’s method to Eq. (22) yields iterations in which the 
following linear systems of equations must be solved: 


VF,, a, s®, am. uPA =-F,, (x, s®, a u®) . (23) 


Henceforth we omit the iteration index k. Assuming that s, u > 0, the Newton 
system (23) can be transformed into a reduced, symmetric saddle point form 


V2,L+ (vn EVA) (Vg) \(Ax\__ [rx (24) 
Vg 0 AJ ha 
i ee u — un 
=: A(x,s,A, U) =!AR =:b 
with diagonal matrix X = diag (4, er Er) and 


ry = VyL + (Vh(x))” ZT! (ye + diag(u)h(x)) , 
ry = g(x). 


The KKT matrices A, arising in every PDIPM iteration k, may become more 
and more ill-conditioned when a KKT point is approached. This is caused by 
the sequence of scaling matrices ©“). If strict complementary holds at the KKT 
point x*, if x — x* and if h;(x*) is an active inequality, one can show that 
5 > oo. On the other hand, 5“ — 0, if h;(x*) is inactive [14].? For this 
reason, many IPM software packages (like IPOPT [7]) use direct methods such as 
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LDL! -factorizations to solve the appearing linear systems. These methods are less 
sensitive to ill-conditioned matrices. 

In contrast, iterative linear solvers like GMRES [15] are very sensitive to ill- 
conditioned matrices. A good preconditioner is needed to lower the matrices’ 
condition numbers. In exchange, they offer a higher potential of parallelization and 
allow the use of inexact Interior Point methods, see Sect. 3.3. 

The distribution of the spectrum of a matrix K, defined as 


o(K):= {A € C, Ais an eigenvalue of K}, 


is an indicator of the convergence behaviour of GMRES applied to Kv = b. By 
rule of thumb, a spectrum which is clustered away from 0 is beneficial to the rate of 
convergence of GMRES. [11, Th. 4.62] 

For general non-linear optimisation problems with the corresponding KKT 
matrix 


T 4 x p. 
c= (PR), mente, nem mem, (25) 


the so-called constraint preconditioner is an example of a preconditioner. It is given 


by 
P H BT 
= (26) 


with H being an approximation to H, that is easy to factorize, e.g., H= diag(H) 
[14]. 

One can show that o(K~!K) = {1} U õ with |ö| = n* — nè. Therefore, at most 
n* —n* eigenvalues of K!K are not equal to 1. The distribution of these remaining 
eigenvalues depends on how well H approximates H on the nullspace of B [10].° 

In Sect.4, we describe how to take advantage of the special structure given 
by (16) to construct a parallel preconditioner. Furthermore, we will see that 
there is a close relation between our proposed preconditioner and the constraint 
preconditioner defined above, see Sect. 4.4. 


Proof (sketch) The KKT conditions require that uf (h;(&*) + s*) = 0. If hj(x*) is active, then 
hj(x*) = 0 and s* = 0. Upon strict complementarity follows that u} # 0. Finally, if x) > x*, 


then uP ur, a — s} and yO = b 


3Let Z denote the nullspace of B, ie. BZ = 0. Then H approximates B on its nullspace, if 
(ZTAZY KZTHZ)*1. 
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3.3 Inexact PDIPM 


In contrast to direct methods, iterative solvers allow solving linear systems with 
prescribed accuracy. One can exploit this fact by means of inexact Interior Point 
methods. Within these methods, the linear systems corresponding to the first PDIPM 
iterations are solved with a rather low accuracy and therefore a low number of 
iterations. As the PDIPM iterate approaches the exact solution, the accuracy of the 
linear solver is increased. We use the approach from [3] for determining the accuracy 
in each PDIPM iteration. Here, (23) has to be solved with residual F (ky 1.6; 


VE), x, 5,4, DAV = Fa, 5, 10, pV) 
+7, 
a x a (KT 5 (k) : R a š . 
that satisfies Ir Io < ne with forcing sequence 7x € (0,1). With this 
choice, the update step A® satisfies 


VFo(x™, s®, a”, py A SENI Foe, s®, ICH uw) 


+ r® 


with |r®]a < (ne + MW ||Foa®,s®,ı®, w)||2. Therefore, “yk + ñk can be 


viewed as forcing sequence of inexact Newton methods” [3]. In our numerical 
experiments, we set 7x = 0.1 for all k and computed the relative tolerance for 
the iterative solver as follows 


pr sb) 


ee 0.1 —— m. 
nt Fy, 7, s®, A0, p®)|| 


We stopped to decrease the relative tolerance when it reached a value smaller than 
10710, 


3.4 Termination Criteria for PDIPM 


We use the following criteria for monitoring the progress of PDIPM [21]: 


1 + max{||x |l; lis llo} 


c (x sS À u) Pan |V. L(x, A, 1) loo (28) 
gradio PAE TE max{llàloo, lalo} 


C feas (X, S, à, u) := (27) 
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T 
Ccomp(X, 5, A, 4) := _ (29) 
re ao (30) 


with x” denoting the previous iterate. 


4 Schwarz Domain Decomposition Method for DOPF 


The key idea behind Schwarz domain decomposition methods applied to DOPF 
problems is the decomposition of the set of time steps T into several smaller subsets. 
Due to this decomposition, it is possible to define corresponding submatrices of 
the global KKT matrix defined in (24) that are smaller in size and therefore faster 
to factorize. Furthermore, the factorization of these submatrices can be done in 
parallel. After computing subsolutions corresponding to these submatrices, one can 
reconstruct an approximate solution to the global system (24), which is used as 
preconditioner within a Krylov-type iterative solver. 

In Sect. 4.1 we introduce the mathematics, which describe the additive Schwarz 
Method (ASM) in the context of DOPF problems. In Sect. 4.2 we briefly describe 
some issues related to its implementation. Section 4.3 is intended to provide a 
deeper insight into the subproblems that have to be solved when the ASM is 
applied. Finally, we describe the relation between the ASM and the above mentioned 
constraint preconditioner in Sect. 4.4. 


4.1 Mathematical Formulation of the Additive Schwarz Method 


For the application of the ASM as a preconditioner for the KKT matrix 
A(x, s, Au) € ROX) defined by (24), it is needed to decompose 
the set of time steps T = {1,..., Nr} into q non-overlapping subdomains: 


with T; NT; = Ø fork #1 and Ty := {7,7 +1,...,é°}. 


Afterward, each subdomain T; is expanded by additional s,; time steps on both ends, 
yielding an overlapping decomposition of T (see Fig. 1): 


q 
T(J TE ee 
I=1 
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Ti-1 Tizi 
pom > 
-> 0-0-0-0-00-0-00-00-0-00-6 0-0- 
LS 
Tı 


Fig. 1 Decomposition of time steps with overlap so; = 1 


with 


T — so, I>1 P tf +8, l<q 


E, l=1 ae l=q 


Typically, so; € {1, 2, 3}. 
Furthermore, we use the definitions of Sect. 2.1.3 to introduce 


at 
teT) 
nı := n} + u; 
n:=n"+n*. 
Henceforward, all expressions involving” refer to the non-overlapping subdomains 
T;. 
When restricting the optimisation variables to the components, which belong to 
T;, we write 


x = [x ]rer, and (31) 


rt 
am = [Y ier, = ) . (32) 
N 


Ax Ax] 5 
A = = e R”. 33 
RW (a) 7 ae (33) 


This vector contains the components of the solution vector of (24) belonging to the 
time steps in T;. Moreover, let R; € {0, 1}"*” be a restriction matrix corresponding 
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to the subset T; defined by its action: 


Ax Ax R* 0 
R = , where R; = l : 34 
i (2) (2) z] nn | 0 a) = 


The matrices Rj and Ri are composed of either zero or identity blocks. As an 
example, consider the case of! Æ 1 and/ Æ q: 


Ri = (0 In; 0) € 10,7", 
_ (35) 
RÈ = (0 1,30) € 0,1". 


Note that the size of the zero matrices varies with each restriction matrix and 
each /. Analogously, we define R; as restriction operator corresponding to the non- 
overlapping subdomain T. 

With these definitions at hand, it is possible to extract submatrices A; from the 
KKT matrix A by multiplying it with Rz: 


Aj i= RıAR]. (36) 
In Sect. 4.3, we explore the structure of these submatrices. 
On the important assumption, that all submatrices A; are non-singular, we are 


able to formulate the multiplicative Schwarz Method (MSM) for approximately 
solving AAr = b. We do this in the form of an algorithm [18]: 


Algorithm 1 Multiplicative Schwarz method 
for /=1,..., q do 
I 1-1 Ty-1 i-1 1-1 
AL = AL! + RTA Rı(b- AAR )= 45 +ô 
end for 
q 
Return Ak 
In every iteration /, the current error with respect to the exact solution Ar, given by 
el = Ar- Aj’, 
satisfies 


Ae’! = r!=! with residual vector r7! = b — AAR! ` 
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By restricting r/~! to subdomain T; and solving with Ay’, one obtains an 
approximation 


ei := A7 Rır'! 


to the exact error Rje'=! on subdomain T}. Prolongation with RT yields a correction 


ö = RF é! to the current approximation anf Thus, the multiplicative Schwarz 
method can be seen as “defect correction algorithm”.* 
Unfortunately, the the MSM is a sequential algorithm. In order to parallelize it, 


one omits residual updates in each iteration, yielding the ASM [18]: 


Algorithm 2 Additive Schwarz method A4, = ASM(b) 


Set 49 =0 
for Tel... do j 
1 _ ,l- T ar 
Al, = AR! + RTA! Rib 
end for 
Return A% 


which can be written as 


q 
A = Masmb with Masy := x RFAT'Rı. 
I=1 


The right preconditioned linear system for solving AAr = b is then given by 
AMasmu=b, Ar = Masmu. (37) 


Since Masm is symmetric but in general not positive definite, we use the Gen- 
eralized Minimal Residual (GMRES) method for solving the non-symmetric sys- 
tem (37). In general, the ASM computes a less accurate approximation to the 
solution of AAr = b and therefore needs an increased number of GMRES iterations 
compared to the MSM given by Algorithm 1. However, the ASM offers a higher 
potential of parallelization, which results usually in better parallel scaling. 


*In contrast to “classical defect correction” algorithms, there is no converging sequence 
(AR), = — AR. Nonetheless, the algorithm’s construction implies that Ad, SAF, 
= 
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4.2 Implementation 


In our implementation, we use one processor per subdomain T; and distribute A 
such that every processor stores RA in its local memory. R/A contains the non- 
overlapping portion of rows of A; = RıART ‘ 

To set up Masm once, every processor first has to form its local submatrix Az, 
i.e., in this step the overlapping part of A; has to be communicated to processor / by 
processor l — 1 and/ + 1. Afterward, each processor computes an LU -factorization 
of A), i.e., Ay = LıU). This step doesn’t involve any inter-processor communication 
and can be done in parallel. 

The employment of the ASM as a preconditioner inside the GMRES method, 
requires to compute M4s5,yv“ in each iteration k. v“ denotes the k-th basis vector 
of the Krylov subspace .%,(AM, sy, b). On the assumption that the basis vector’s 


data layout is the same as the matrix’s one, i.e. every process stores a? = Rv, 

communication with process / — 1 and / + 1 is required to multiply Masm with 

v® > The computation of Am is done by one forward substitution with L; 

and a backward substitution with U, respectively. This step does not involve any 

1, (k) 
l 


communication. After this, the local solution A, v; is prolonged back to obtain 


the global result of the matrix vector product MA suv, This again requires 
communicating the overlapping part of uf to process l — 1 and/ + 1. Finally, it’s 
necessary to multiply A with Masu v Due to A’s data layout and its off-diagonal 
elements, corresponding to intertemporal couplings, additional communication is 
required. 

One can further improve the performance of ASM by using the so-called 
restricted version of ASM [4], which is given by 


P 

M; ASM = So RIARI. (38) 

I=1 
For this preconditioner, just the non-overlapping part of the local solution um is 
prolonged instead of the entire (overlapping) vector. Experiments show a beneficial 
behaviour in terms of GMRES iterations compared to standard ASM [4]. Further- 

more, prolongation by Rı doesn’t involve any communication. 

To further stabilise the LU-factorization, a small regularisation parameter e = 

1078 is added to the KKT matrix A, i.e. A + el. 


5The reader may wonder why we have to make an assumption about the basis vectors’ data layout: 
Our GMRES implementation is part of the external library PETSc and a quick glance at their 
source code did not reveal to us how they manage the according data. 
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Fig. 2 Permuted KKT matrix A with Nr = 12, sy) = 1, q = 4 and A; = PA; PT 


4.3 Structure of Local KKT Matrices 


In this section, we investigate the structure of the local submatrices A; defined in 
Eq. (36). A good way to get an idea of their structure is to permute the KKT matrix 
A in such a way that all rows and columns belonging to one time step are grouped 
together. This yields a block diagonal matrix with some off-diagonal elements. 
Every block represents one time step and the off-diagonal elements arise due to the 
time-dependent constraints. Extracting the time steps belonging to the subdomain 
T; results in a permuted version of the submatrix A7.° Its structure is exactly the 
same as the one of the original KKT matrix A. If the off-diagonal elements are 
taken into account as “boundary conditions”, it is possible to interpret the submatrix 
A, as representing an optimisation problem in its own, i.e. the original dynamical 
OPF problem is reduced to the subdomain T;. Figure 2 tries to demonstrate this. As 
a consequence the submatrices A; tend to be ill-conditioned. 


Note that this extraction cannot be done by multiplying A with R; and RT. 
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4.4 Relationship Between the ASM and the Constraint 
Preconditioner 


To see a relation between the additive Schwarz method and the constraint precon- 
ditioner, as introduced in Sect.3.2, it is necessary to consider a dynamical OPF 
problem without time-dependent equality constraints. In our case, we had to drop 
the storage facilities 85 to do this. Furthermore, it is important to consider a non- 
overlapping decomposition of the time domain, i.e. So = 0. In this case, the 
corresponding ASM preconditioner, 


q 
MASM = > RT A 'Rı, (where Al = RıART e RU*") (39) 
I=1 


reduces to a block Jacobi method with omitted couplings between variables 
belonging to a and tay for each subdomain /. Since there are only time-dependent 
inequality constraints left, the coupling between the variables belonging to different 
time steps has to arise in the (1, 1)-block of the KKT matrix, namely in the 
(Vh)! 5(Vh) part. This means that multiplying Masm with b is like solving a 
modified KKT system, or equivalently, 


` -1 
7 T 
MASM = G i ) ; (40) 


H denotes the KKT matrix’s modified (1, 1)-block. Thus, M asm has the form of a 
constraint preconditioner. 

At this point, it also becomes clear that the ASM can only be interpreted as 
a constraint preconditioner, if there are no time-dependent equality constraints. If 
there existed time-dependent equality constraints, the (2, 1)- and (1, 2)-block of 
Mis m would have changed too. 

As pointed out in Sect. 3.2, 1 is an eigenvalue of MasmA of multiplicity 2n* and 
the remaining n* —n? eigenvalues are solutions of a generalized eigenvalue problem 
of the following form [10]: 


ZT (v2 L+ (Wh)? E(VhA))Zv = AZ AZv. (41) 
nd 


=H 


For a low number of subdomains q, one can expect H to bea good approximation 
to H, leading to a clustered spectrum o(M AsmMA) close to one. However, the 
more subdomains, the more neglected couplings and the less accurate does H 
approximate H. This results in a more scattered eigenvalue distribution of Mas A, 
but still a large number of eigenvalues are equal to 1. 
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5 Numerical Experiments 


In this section, we present some results for our proposed method applied to aDOPF 
problem. In its first part, we discuss the parallel speedup of our solver. In the second 
one, we investigate the way in which the KKT matrix’s spectrum o (A) changes 
with the iterations of the interior point method. The test grid, which we used for our 
experiments, represents a possible variant of the German Transmission Grid in the 
year 2023 and consists of 1215 nodes, 2247 AC lines, 8 DC lines, 181 conventional 
power plants, 355 renewable energy sources and 67 storage facilities. For more 
details, we refer the reader to [12]. 

For solving (16), we use the PDIPM algorithm mips which is written in MATLAB 
code and part of MATPOWER [21]. In this algorithm, we replace the standard 
MATLAB backslash operator \ for solving linear systems by our own linear solver. 
Our solver applies a right preconditioned GMRES method. We use the ASM as 
the preconditioner. For computing the LU-factorizations of the local systems Aj, 
we use the software package Intel® MKL PARDISO. Our solver is written in C++ 
and makes use of the KSPGMRES and PCASM methods provided by PETSc [2], 
which is compiled in Release mode with the Intel compiler 17.0. The computation 
of the eigenvalues of the KKT matrices and the preconditioned KKT matrices have 
been done with SLEPc [9] and MATLAB. Inter-process communication is realised 
by means of the Message Passing Interface (MPI). The numerical experiments were 
carried out on the BwForCluster MLS&WISO Production at Heidelberg University 
with two Intel Xeon E5-2630v3 processors i.e. 16 CPU cores at 2.4 GHz per node, 
and 64 GB working memory. 


5.1 Parallel Speedup 


In our numerical experiments, we concatenated the 48 h load profile, which was part 
of the described test grid, to obtain a test case comprising 48 days with a temporal 
resolution of Ih, i.e. Nr = 960 time steps. This causes the linear system (24), 
which has to be solved in every PDIPM iteration, to have the dimension n* + n? ~ 
6.168 - 10°. 

Furthermore, we set a = 0.4, which led to ramping constraints limiting 
the change of power generation of conventional power plants by 40%/h of their 
respective maximum power generation. 

We set €feas = €grad = €cost = €comp = 1075 as PDIPM termination 
criteria and limited the number of PDIPM iterations to 300. The number of GMRES 
iterations was restricted to 1500. The GMRES method was allowed to restart after 
300 iterations. Its relative residual tolerance rtol was determined as described in 
Sect. 3.3. 

The ASM method was not restricted and the overlap so; was set to 3 for all tests. 
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Fig. 3 Speedup up to the 100 
iteration on the x-axis for —— ASM: #2 
different numbers of —— ASM: #16 
processors N ASM: #64 
= —— ASM: #192 
a c | — ASM: #320 
3 
3 10 ----LU: #16 
a. 
a Na 
2 nr ee 
IS = 
10 20 30 40 
PDIPM iteration k 


The parallel speedup on g processors for the first k PDIPM iterations is defined 
as sq(k) := TE We computed the sequential reference time Tı(k) by solving 
the linear systems (24) using a LU -factorization and adding up the times needed 
to solve them up to the k-th PDIPM iteration. MKL PARDISO’s sequential LU- 
factorization algorithm was applied. 

Analogously, 7,(k) was obtained by adding up the times needed to solve the 
same linear systems, but this time iteratively, i.e. the GMRES method and the 
ASM were used. In Fig.3 we plotted the results of our speedup computations for 
different numbers of processors. The dashed graph represents the speedup obtained 
when the KKT systems (24) are solved directly with MKL PARDISO’s parallel 
implementation of the LU-factorization using 16 processors. For this amount of 
processors we observed the best speedup. 

Besides, we scaled the y-axis logarithmically to ease drawing a comparison with 
the LU-solver and to emphasise the large speedups during the first iterations.’ We 
observed for our proposed iterative solver, no matter the number of processors 
q, a clear speedup in comparison with the LU-solver. The parallelization, due to 
the ASM, shows its power during the first PDIPM iterations. Unfortunately, this 
initial speedup decreases rapidly. After 40 iterations we end up with a total speedup 
of five. This decrease is accompanied by an increase in GMRES iterations. To 
demonstrate this, we plotted in Fig. 4 the number of GMRES iterations needed to 
solve the KKT system at the k-th PDIPM iteration. The first thing to note is, that 
the number of GMRES iterations does not only change with PDIPM iterations, but 
also with the number of processors used to solve the KKT systems. It’s likely that 
the increase in GMRES iterations with PDIPM iterations is caused by a change 
in the spectrum of the original KKT matrix A. The increase with the number of 


7We remark that the PDIPM needed in most of the tested cases more than 40 iterations to converge, 
yet, it converged after at most 45 iterations. 
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Fig. 4 The number of 700 
GMRES iterations needed to 2 
solve the KKT system © 600 
corresponding to the PDIPM 5 500 
iteration on the x-axis 2 400 
i 
& 300 
5 200 
+ 100 
10 20 30 40 
PDIPM iteration k 


processors can be explained as follows: The more proccesors we use, the smaller 
the subproblems, described by A;, become and, hence, the more intertemporal 
couplings are neglected, which reduces the “approximity” of Masy to A~!. And 
the closer Masm is to A~!, the more the spectrum of the preconditioned operator 
o(AMasm) is clustered around 1. 

Before we start to actual investigate the spectrum of the matrices in the next part 
of this section, we would like to summarise our findings. Even though we observed 
a rapid decrease in the speedup, it was possible to obtain a moderate speedup with 
only 16 processors and in comparison to a standard parallel LU-factorization also 
using 16 processors, our solver was more than a factor 2 faster. 


5.2 Eigenvalue Distribution 


Since the GMRES method converges faster if the eigenvalues of matrix A are 
clustered around 1, i.e. it needs less iterations, we would like to investigate how 
much the ASM, applied as preconditioner, improves the eigenvalues’ distribution. 
Therefore, we calculated the eigenvalues with the largest and smallest magnitude of 
the KKT matrix A and the preconditioned KKT matrix AMAsm, respectively. We 
write A), = max{|A|: A € o(A)} and A/,, = min{|A|: A € o(A)}. Seeing that it 
takes very long to compute, or that it is even not possible to compute the eigenvalues 
for the complete Nr = 960 time steps, we decided to calculate them for Nr = 96 
time steps and for each PDIPM iteration. We set the overlap to sy; = 1 and used 
q = 32 subdomains. The results are presented in Fig.5. The red lines depict the 
development of yn ax and 24 in Of the KKT matrix, respectively. Thus, they form 
the limits for all possible magnitudes of eigenvalues corresponding to A. The latter 
is indicated by the red shaded area. Analogously, the blue lines and the blue area 
represent the spectrum of the preconditioned KKT matrix. 

The spectrum of the preconditioned KKT matrix, is clustered around 1 in every 


interior point method iteration. Looking at the magnified part of the graph shows 
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Fig. 5 Variation of the 
spectrum of the original KKT 
matrix and the preconditioned 
one 


An in/max 


# PDIPM iterations 


that en becomes smaller between the 10-th and the 20-th iteration. This may 
explain the increase in GMRES iterations observed in our numerical experiments 
and plotted in Fig. 4. 


6 Conclusion 


In this work we proposed a way of solving linear systems arising from Dynamic 
Optimal Power Flow (DOPF) problems in parallel by means of an overlapping 
Schwarz domain decomposition method, namely the additive Schwarz method. 
It was shown how to apply this method in the context of DOPF problems and 
that the obtained submatrices correspond to localised formulations of a DOPF 
problem with additional boundary conditions. Numerical tests on one large-scale 
DOPF problem showed that the combination of the Generalized Minimal Residual 
(GMRES) method and the Additive Schwarz Method (ASM) is able to solve DOPF 
problems. Furthermore, we were able to show that this combination yields good 
parallel speedup for some of the PDIPM iterations. And, in a small example, we 
demonstrated that the ASM is a good preconditioner in the sense, that it clusters the 
eigenvalues around 1. 

Unfortunately, the proposed solver had problems with some of the KKT systems 
arising in the large-scale DOPF problem. Therefore, it is necessary to improve it 
and this will be one of our main goals in future work. 
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Part II 
Energy System Integration 


Optimal Control of Compressor Stations A 
in a Coupled Gas-to-Power Network gai 


Eike Fokken, Simone Göttlich, and Oliver Kolb 


Abstract We introduce a tool for simulation and optimization of gas pipeline 
networks coupled to power grids by gas-to-power plants. The model under consid- 
eration consists of the isentropic Euler equations to describe the gas flow coupled 
to the AC powerflow equations. A compressor station is installed to control the 
gas pressure such that certain bounds are satisfied. A numerical case study is 
presented that showcases effects of fast changes in power demand on gas pipelines 
and necessary operator actions. 


Keywords Coupling of gas and power networks - Compressor stations - Optimal 
control 
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1 Introduction 


Renewable power sources have an ever increasing share of all power sources. 
Though renewable energy has been developed in recent years with great success, 
its intermittent and unpredictable nature raises the difficulty to balance the energy 
production and consumption [6, 18]. A frequent proposal is to use gas turbine plants 
to compensate for sudden drops in power of renewable sources because these plants 
are relatively flexible in comparison to coal or nuclear plants. A welcome advantage 
of this approach is the possibility to run gas plants with fuel produced from 
renewable electricity via power-to-gas plants, thereby reducing or even negating 
carbon emissions of the gas plants. For a review of power-to-gas capability see [6]. 
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It is desirable to have a joint optimal control framework for power and gas sector 
of the energy system to model this compensation. So far only steady-state flow in 
the gas network has been considered [3, 18, 20], which may be too coarse for several 
applications. Therefore, we focus on an optimal control strategy for the instationary 
gas network model [1] coupled to a power grid [2] via compressor stations, see for 
example [8, 15, 17]. The mathematical foundation for the gas-to-power coupling 
has been recently introduced in [5], where conditions for the well-posedness have 
been derived and proved. This work is the first attempt to model this interaction 
and yields an understanding of the underlying equations. The next step will be real- 
world scenarios. 


2 Optimal Control Problem 


The gas dynamics within each pipeline of the considered gas networks are mod- 
eled by the isentropic Euler equations, supplemented with suitable coupling and 
boundary conditions. For the power grid, we apply the well-known powerflow 
equations. The coupling between gas and power networks at gas-driven power 
plants is modeled by (algebraic) demand-dependent gas consumption terms. To 
react on the demand-dependent influences on the gas network, controllable devices 
as compressor stations are considered within the gas network. The aim is to fulfill 
given state restrictions like pressure bounds whereas at the same time the entire fuel 
gas or power consumption of the compressor stations is to be minimized. 

The network under consideration is similar to the one presented in [5] and is 
depicted in Fig. 1. In addition, there is now a compressor with a time-dependent 
control u(t), which is used to satisfy pressure bounds. The optimal control problem 
is to minimize compressor costs while satisfying power demand and gas dynamics: 


min,(r) Compressor costs 
s.t. isentropic Euler equations (Sect. 2.1) 
gas coupling conditions (Sect. 2.2) 
compressor equation (Sect. 2.3) 
powerflow equations (Sect. 2.4) 
gas-power-coupling (Sect. 2.5) 
pressure bounds. 


Mathematically, this is an instationary nonlinear optimization problem con- 
strained by partial differential equations, see [9] for an overview. To solve the 
problem, we make use of a first-discretize-then-optimize approach and apply the 
interior point solver IPOPT [16]. The necessary gradient information for IPOPT, i.e., 
gradients with respect to all controllable devices, is efficiently computed via adjoint 
equations. Here, the underlying systems can be solved time-step-wise (backwards 
in time), where additionally the sparsity structure is exploited. 
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Fig. 1 Gas network connected to a power grid. Red nodes are PQ/demand nodes, green nodes are 
generators (PV nodes) and the blue node is the slack bus (also a generator, with gas consumption 
of the form e(P) = ag +a; P + a P?). The circle symbol indicates a compressor station 


We remark that the cost function is given in Sect.2.3, while the bounds on 
the pressure are introduced as box constraints within the numerical optimization 
procedure in Sect. 3. Within the Sects. 2.1, 2.2, 2.3, 2.4, and 2.5, we now describe 
the constraints of the optimal control problem in detail and focus on the technical 
details. 


2.1 The Isentropic Euler Equations 


The gas network is modeled by the isentropic Euler equations (see [1, 4]), which 
govern gas flow in each pipeline between nodes, 


p q 0 
+ = 5 1 
a la + 2) ee n m 


where p is the density, q = pv is the flow, p(p) = xp” is the pressure function. In 
our example we use y = | and «k = e = (3405), that is, the isothermal Euler 
equations with speed of sound c. The Euler equations must be met for 0 < x < le 
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and t > 0, where le is the length of the pipe. Furthermore, S is a source term, 


Ag)alal 
= 
e 


where A is the solution to the Prandtl-Colebrook formula, 


a 251, ke 


The Reynolds number Re is given by 


d 
Re(q) = =q 
n 
with dynamic viscosity 
k 
ei 
ms 


The roughness ke and diameter de of the pipes are all the same in our example and 
given by 


ke = 5-107*m, 
de = 6< 10" m. 


2.2 Coupling at Gas Nodes 


At the nodes we use the usual Kirchhoff-type coupling conditions: The pressure is 
the same near the node in all pipes connected to it and the flows must add up to zero 
(where the sign for inflow is positive, the sign for outflow negative), 


Ppipe = Pnode» (2a) 
> Ipipe = > Ipipe- (2b) 
ingoing pipes outgoing pipes 


The example gas network we are using is a small part of the GasLib-40 
network [10]. In Table 1 the only remaining parameter, the length le of each pipe, 
is gathered. Note that no length is given for the arc connecting SO and S17, because 
only a compressor is situated between these nodes. 
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Table 1 Parameters of the Pipe [From [to | Length Ze [km] 
pas network P10 |S4 |s20 | 20.322 
P20 |ss Iso | 20.635 
P21 |sır |s4 | 10.586 
p22 |sır |ss | 10.452 — 
P24 |Ss8 | $20 | 19.303 
P25 |s20 | $25 | 66.037 


2.3 Compressor Stations 


To compensate for pressure losses in the gas network, we consider compressor 
stations, which are also modelled as arcs. Those arcs have (time-dependent) in- 
and outgoing pressure (Pin, Pout) and flux values (gin, Gout). In general, the two 
separate flux values allow the modelling of fuel gas consumption of the compressor 
station, whereas we will consider an external power supply for the compressor and 
therefore have gin = Gout- The power consumption is modelled as a quadratic 
function of the power required for the compression process. Denoting this as 
function P(pin(t), gin(t), Pout(t), Gour(t)), our objective function in the optimal 
control problem is of the form 


T 


J P (pin(t), qin(t), Pout(t), Jout(t))dt. (3) 


0 


Note that the power consumption does not influence the network dynamics and 
is therefore only of interest for the optimization procedure. For our investigations 
below it is sufficient to know that the consumption and therewith the costs increase 
if the ratio Pout/ Pin increases. For the details of the power consumption model, we 
refer to [17, Section 3.2.3]. The influence of the compressor station on the network 
dynamics is modelled by the control u(t) of the pressure difference: 


Pout(t) — Pin(t) = u(t). 
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2.4 Power Model 


For the power grid we use the AC powerflow equations (see [7] for an introduction), 


N 
Pk =) \|Vel|V;| (Ge; cos(Ør, j) + Bey sin(@x j)), 
j=l 


(4) 
N 
Ox = ) \Vel|V;| (Gaz sin(bx,j) — Bej cos(x,)), 


k=1 


where Px, Ox are real and reactive power at node k,|V;| is the voltage amplitude, Øx 
is the phase (and x, ; = k — $j) and Byj, Gx; are parameters of the transmission 
lines between nodes k and j or of the node k for Bg and Gxx. Each node is either 
the slack bus (in our case N1; V and & given), a generator bus (N2 and N3; V and 
P given) or a load bus (N4 through N9; P and Q given). All in all for N nodes we 
get 2N equations for 2N variables. 

The considered power grid is taken from the example “case9” of the MAT- 
POWER Matlab programming suite [19]. A per-unit system is used, whose base 
power and voltage are 100 MW and 345kV respectively. The corresponding node 
and transmission line parameters are gathered in Table 2, these are the entries of the 
nodal admittance matrix (see [7]). 


2.5 Coupling 


The last ingredient is a model for converting gas to power at a gas power plant. In 
our example, it will be situated between the nodes S4 and N1 and convert a gas flow 
£ into a real power output P according to 


e(P) = ao +a P + aP’. (5) 


The flow £ must be diverted from the pipeline network, hence the coupling condition 
at node S4 must be changed to 


Pin = Pout, 
(6) 
qin = Jot + € . 


The details of simulation of such a combined network and the treatment of all arising 
mathematical issues can be found in [5]. We now showcase a concrete example. 
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Table 2 Parameters of the 
power grid (p.u.) 


3 Numerical Results 


3.1 Problem Setup 


(a) Busses 


Node 
N1 
N2 
N3 
N4 
N5 
N6 
N7 
N8 
N9 


(b) Transmission lines 


Edge 
TL14 
TL45 
TL56 
TL36 
TL67 
TL78 
TL82 
TL89 
TL94 


G 

0.0000 
0.0000 
0.0000 
3.3074 
3.2242 
2.4371 
2.7722 
2.8047 
2.5528 


From | To 


N1 
N4 
N5 
N3 
N6 
N7 
N8 
N8 


[N9 


N4 
N5 
N6 
N6 
N7 
N8 
N2 
N9 
N4 


B 


73 


—17.3611 
—16.0000 
—17.0648 
—39.3089 
—15.8409 
—32.1539 
—23.3032 
—35.4456 
—17.3382 


G 
0.0000 
—1.9422 
—1.2820 
0.0000 
—1.1551 
—1.6171 
0.0000 
—1.1876 
—1.3652 


B 
17.3611 
10.5107 
5.5882 
17.0648 
9.7843 
13.6980 
16.0000 
5.9751 
11.6041 


As already noted, we consider a small part of the GasLib-40 network from [2] 
consisting of 7 pipelines with a total length of 152km. This network is extended 
by a compressor station and additionally connected to a power grid with 9 nodes 
by a gas-to-power generator. For this coupled gas-power network, we simulate a 
sudden increase in power demand within the power grid and study its effect on the 
gas network. The considered compressor station is supposed to compensate part of 
the pressure losses in the gas network such that a given pressure bound is satisfied 
all the time, while power consumption of the compressor is minimized. 
To complete the problem description, the following initial and boundary condi- 
tions are given: 


P(t) and Q(t) at load nodes, 


P(t) and V(t) at generator nodes, 


V(t) and @(f) at the slack bus, 
p(t) at S5, 

q(t) at S25, 

P(x, 0), q(x, 0) for all pipelines. 
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Table 3 Initial conditions of 
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Fig. 2 Power and reactive power at demand node N5 (above) and the slack bus (below) 


More precisely, the initial conditions for the power network are given in Table 3. 
These remain constant over time except for the power demand at node N5, which 
changes linearly between t = 1h andt = 1.5h from 0.9 p.u. to 1.8 p.u. for the real 
power and from 0.3 p.u. to 0.6 p.u. for the reactive power, see also Fig. 2. 

For the gas network the incoming pressure at S5 is fixed at 60bar, the outflow 
at S25 is fixed at q = 1007 . RX, where 00 = 0.7855 and Ae = rh. The 
fuel consumption parameters we use in equation (5) are given by ap = 2, aj = 5, 
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az = 10. Since the data for the considered gas and power networks are taken from 
different sources, the parameters of the gas-to-power generator are chosen in such a 
way that a significant influence is caused. 

Further, the pressure at S25 is supposed to satisfy a lower pressure bound of 
4 bar, i.e., 


Ps25(t) > 41bar 


for all times t, where we consider a time horizon of T = 12h. 

As we will see below, the compressor station between nodes SO and S17 will 
have to run at a certain time, i.e., u(t) > 0, to keep this pressure bound. In general, 
a solution to the described optimal control problem consists of the control u(t) and 
the entire network state for all times f, i.e., 


— P(t), Q(t), V(t), &(t) for all nodes in the power grid, 
— p(x,t), g(x, t) for all pipelines in the gas network, 


fulfilling the given model equations and the pressure constraint. 


3.2 Discretization and Optimization Schemes 


To solve the described optimal control problem, we follow a first-discretize- 
then-optimize approach. The model equations of the power grid only require a 
discretization in time, which means that the given boundary conditions and the 
powerflow equations (4) hold for discrete times t; = j At with At = 15min and 
j € {0, ..., 48} in our scenario. The discretization of the isentropic Euler equations 
within the pipelines of the gas network additionally requires a spatial grid (here 
with grid sizes Ax, ~ 1 km) and an appropriate discretization scheme. Here we 
apply an implicit box scheme [14], which allows the application of large time steps 
as At = 15min for the considered spatial grid. Considering the isentropic Euler 
equations as a system of balance laws of the form 


ye + fO) = 8), 


the applied scheme reads 


ro 4 rn Y? 2h. y? 


2 2 
(sort) - fo) 


he +87) 
m. 
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The numerical approximation of the balance law is thought in the following sense: 


Y? ~yQ,t) for xe [G- 5) Ax, (J + Ņ)4x), 
t € [nAt, (n+ 1)Ar). 


Together with the algebraic equations modelling the compressor station and the 
coupling and boundary conditions, the discretization process results in a system of 
nonlinear equations for all state variables of the coupled gas-power network. For 
simulation purposes, the entire discretized system is solved with Newton’s method. 
Note that the system can be solved time-step per time-step and that we exploit the 
sparsity structure of the underlying Jacobian matrices. 

So far we can only compute the state of the considered gas-power network (for 
a given compressor control u(t)) and evaluate quantities of interest like the power 
consumption of the compressor station or the pressure constraint within the time 
horizon of the simulation. For the given time discretization, the compressor costs (3) 
are approximated by the trapezoidal rule and formally contained in J(u, y(u)) 
below. Next, we want to solve the (discretized) optimal control problem, i.e., to 
find control values u(t;) such that the pressure constraint is satisfied, while the 
power consumption of the compressor is minimized. For this purpose, we apply 
the interior point optimization code IPOPT [16] to the (reduced) discretized optimal 
control problem. In addition to our simulation procedure to evaluate quantities of 
interest for a given control, IPOPT further requires gradient information of those 
quantities with respect to the control. Based on the considered discretization, such 
information can be efficiently computed by an adjoint approach, which we briefly 
describe in the following. Therefore, we consider the discretized optimal control 
problem in the following abstract form: 


min J(u, y(u)) 

st. ue R", y e RY, 
model equations: E(u, y(u)) = 0, 
constraints: g(u, y) > 0. 


The vector u contains all control variables (here the compressor control u(t;) for 
all times t; = j At € [0, 24]) and y contains all state variables of the coupled gas- 
power network of all time steps t;. The mapping u —> y(u) is implicitly given by our 
simulation procedure. Thus, IPOPT does not have to care about the model equations 
formally summarized in E(u, y), but only about the further constraints g(u, y) and 
minimizing the objective function J(u, y). Accordingly, we need to provide total 
derivatives of J and g with respect to u. 

In the following, we consider the computation of these derivatives only for the 
objective function, since the procedure is identical for the constraints, and we follow 
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the description given in [12]. First of all, the chain rule yields 
d ð 0 d 
— J(u, yu)) = —J(u, ylu)) + —J(u, yu))— yu). (7) 
du ðu dy du 


While the partial derivatives of J with respect to u and y can be directly computed, 
the derivatives of the states y with respect to the control u are only implicitly given. 
Differentiating the model equations E(u, y(u)) = 0 yields 


i ves ie ees 
— u, yu = u u, yu ay u, yu ae 


and therewith (formally) 


a, ay 8 
Er iz e w 1) a E) (8) 


Even though the partial derivatives on the right-hand-side of (8) can be directly 
computed, one would have to solve Nu systems of linear equations here. Instead of 
that, we insert (8) into (7) and get 


d ð 
a ee) = a, a) 


a a -l a 
= y(u)) (Seu. vw) — E(u, y(u)). 
y dy ou 
—— 
=T 


With the so-called adjoint state & as the solution of the adjoint equation 


a £ ð 2 
(Feu. yw)) = (Zu yw)) 7 (9) 
y dy 


we finally have 


d f) a 
g! ur — J(u, y(u)) +87 —E (u, y(u)). (10) 
u ðu ou 


It is the fact that (9) is a single linear system and has a special structure, which 
can be easily exploited (see for instance [11, 13]), which makes the computation of 
derivatives via the presented adjoint approach very efficient. Nevertheless note that 
for given control variables u one still has to solve the model equations to get y(u), 
before one may compute the gradient information via (9) and (10). 
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3.3 Results 


We first discuss the simulation with inactive compressor. In the course of the 
simulation, due to the increase in power demand at node N5, the power demand 
at the slack bus rises as well and leads to increased fuel demand at node S4. This 
increases the inflow into the gas network, as can be seen in Fig. 3. Also the pressure 
in the final node S25 decreases and violates the pressure bound after approximately 
4h, see Fig. 4. 

After the optimization procedure, the compressor station compensates part of the 
pressure losses in the gas network such that the pressure bound is satisfied all the 
time. Since the power consumption of the compressor station is minimized within 
the optimization, the pressure constraint is active after roughly 4 h (see again Fig. 4), 
i.e., the compressor station applies as little as possible energy. 
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Fig. 3 Inflow at the node S5 
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Fig. 4 Pressure at the node S25 
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Conclusion 


The proposed optimization model allows to predict pressure transgressions within a 
coupled gas-to-power framework. Simulation and optimization tasks are efficiently 
solved by exploiting the underlying nonlinear problem structure while keeping the 
full transient regime. This makes it possible to track bounds much more accurately 
than with a steady state model, thereby achieving lower costs. 
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Utilising Distributed Flexibilities A 
in the European Transmission Grid PEE 


Manuel Ruppert, Viktor Slednev, Rafael Finck, Armin Ardone, 
and Wolf Fichtner 


Abstract In this paper, we investigate the effect of distributed flexibilities on the 
operation of the transmission grid. The flexibilities considered are heat pumps, 
electric vehicles, battery energy storage systems and flexible renewable generation. 
For this purpose, we develop a two-stage approach of first determining an optimal 
electricity market solution considering the optimal dispatch of each generation 
element and flexibility. In the second step we determine the required dispatch 
adjustments due to transmission grid constraints and investigate the effect of 
integrating battery energy storage systems into the adjustable generators to solve 
congestions. In our case study, we investigate the central European transmission grid 
for a scenario based on the Distributed Generation scenario of the Ten-Year Network 
Development Plan for the year 2030. Integrating distributed flexibilities leads to a 
strong increase in the security of supply, while the overall effect on the generation 
adjustment is small. A comparison of the results for an AC and DC formulation 
shows that both approaches differ significantly in individual cases. 


Keywords AC optimal power flow - Battery energy storage systems - DC 
optimal power flow - Distributed flexibilities - Distributed generation - Flexible 
demand - Transmission grid 


1 Introduction 


Today’s change in the electricity system is driven by the decarbonisation initiative 
developed to meet emission reduction targets defined in international agreements 
in order to limit the global temperature increase. To reduce the carbon intensity 
of the electricity generation, an increasing amount of electricity generation from 
renewable sources (RES-E) is being installed throughout Europe. The two most 
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significant sources are wind and solar, both characterised by volatility and their 
spatially distributed generation potential [1]. On the demand side, decarbonisation 
efforts have begun to lead to an increased electrification in the heat and transporta- 
tion sectors. In some countries, rapid market growth of such elements which can 
be categorised as distributed flexibilities can be already be observed today. Future 
energy scenarios hint towards a further increase of these developments until 2030. 
With potentially high concurrencies of individual, decentralised demand patterns, 
uncontrolled or market-driven operation of these flexibilities will increase local as 
well as aggregated consumptions peaks in the power sector. 

These changes pose significant challenges to the operation of the European 
transmission grid. First, balancing supply and demand will require additional 
flexibilities in a system dominated by distributed, intermittent RES-E and new 
sources of electrical demand. Second, RES-E at a specific point in time will be 
largely determined by the predominant meteorological conditions with high spatial 
concentration, increasing average distance between generation and consumption 
and thus increase the stress on the grid infrastructure. Lastly, market requirements 
to create a cost-minimal dispatch might contradict those necessary to enable a 
congestion-free transmission grid. 

In the past, extensive research has been performed on different flexibility types 
concerning the technologies, their modelling and their contribution to renewable 
integration. A comprehensive overview on the current developments in electrical 
energy storage technologies and their potential for application in power system 
operation is given in [2] and [3]. Current reviews of different technological options 
can be found for power-to-heat applications [4], residential photovoltaic-battery 
energy storage (PV-BESS) [5], heat pumps (HP) in smart grids [6], and battery 
electric vehicles (BEV) [7]. Most literature focuses on the local balancing of demand 
and supply, either on a single household level or a community. This neglects the 
implications on the grid level or focuses on the regional effects in micro grids while 
including only a small set of technological options. 

There are numerous models, which take into account flexibilities and simplified 
transmission grid constraints. In most cases, these models have been developed 
to analyse electricity markets, e.g. for investment decision support or operation 
decision support. In most studies liberalised markets in the United States or 
Europe are analysed and grid constraints are represented in a linearised way or 
by import/export constraints between market zones. Hence, most techno-economic 
models only take these constraints into account. 

Few models, which consider flexibilities in a grid context take into account 
the full AC formulation of transmission grid constraints. Besides the group of 
studies and models, which consider a linearised form of grid constraints, some 
models can integrate full AC constraints. However, they are limited to either a very 
short time horizon or their focus lies on transient technical aspects like fault level 
detection, transient stability, harmonic analysis, reliability and power flow. Another 
group of models takes into account many flexibilities and AC restrictions but 
only on a distribution grid level, thus limiting system size and thereby complexity 
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significantly [8]. A comprehensive overview on modelling approaches and their 
consideration of the grid can be found in [9]. 

Another area where transmission constraints are explicitly considered is the 
co-optimisation of transmission and generation capacities. The complexity of the 
resulting problem forces modellers to apply linearised transmission grid constraints 
or reduce the number of modelled flexibilities [10]. Furthermore, such models 
are often only validated on small scale or test grids, making it difficult to assess 
their suitability for large real-world power grids [11]. In most cases, the focus 
lies on investment decisions rather than on utilising operational flexibility in future 
transmission systems. [12] provide an overview of requirements and approaches to 
address this problem. 

In this paper, we present a modelling framework to investigate the effect of 
distributed flexibilities on the transmission grid operation and investigate the benefit 
of utilising distributed flexibilities to reduce grid congestion while taking into 
account the full AC model of the transmission grid. 


2 Modelling the Interaction of Market and Grid Operation 


In order to model today’s interaction between the electricity market and transmis- 
sion grid operation, we developed a two-step approach to model market and grid 
operation. In the first step, we use a linear optimisation approach to determine the 
minimal-cost, copperplate-based results of the interconnected electricity market. In 
the second step, we determine the required dispatch adjustments due to congestion 
in the transmission grid using a multi-period nonlinear optimisation model. In the 
following section, the general modelling of the electricity market is described, 
followed by the integration of distributed flexibilities and the application of the 
method for modelling the transmission grid operation. 


2.1 Electricity Market 


To model the interconnected electricity market consisting of multiple market areas, 
we use a linear optimisation approach with the objective function of minimal 
total annual system cost under the assumptions of perfect foresight and perfect 
competition.! The system cost consist of the aggregated variable cost of electricity 


LA representation of the electricity market using an optimisation approach with cost minimisation 
objective function and perfect foresight does not necessarily allow for a realistic representation of 
market prices, as scarcitiy rents cannot be considered. As the market model is first and foremost 
utilised to obtain dispatch results for the transmission grid simulation in our application, using the 
described approach is a reasonable simplification for this purpose. 
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production Cg for the electricity generation p, of generator g in time step t and 
the cost of load shedding Czs for the amount of load shedding pis of demand d in 
time step t. 


4 LS 
min x Cg * Pg + % CLs * Pat 
" geG,teT deD,teT 


In every market area m, the electricity demand Py; in every time step must be 
balanced out against generation, load shedding (LS) and interconnector flows into 
the market area pe,r. 


XO Pett >) Pat DS Pet = >) PaıYmeM,teT 


geGM deDM ceIc" deDM 


The bounds of each generator are determined by the minimum capacity pee and 
maximum capacity Poe of the respective generating unit. While the boundaries 
of conventional generation are determined by the technical parameters of the 
generator, the boundaries of RES-E are determined by the available generation due 
to the intermittent nature. Available interconnection flows between market areas 
are restricted by the available exchange capacity prin. P:"@*, In order to ensure 
feasibility of the problem, LS up to total demand can be performed at each time 
step. 


pur < Pet < Por YgEG,tET 
PER < pep < POO YceIC,teT 


0< ps < Pas Vd € Dt eT 


2.2 Modelling Flexibilities 


In the following, we apply a set of general equations to model various types 
of time-coupled flexibilities, ranging from different power storage technologies 
to consumption and renewable flexibilities. Ignoring self-discharging losses, the 
energy level vs, of a storage s € S at time step t € T = {1,... NT} equals its 
previous level, the external power inflow qe or outflow cy and the sum of charged 


and discharged power provided by connected generators pg; € GS CG: 


in out 
Us = Usu- D> Per *ng— D> Pens + 04% io 
geGms geGous 


VseS,teT 
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For a simpler notation, the set of storage-connected generators GS is split such 
that G°"'S denotes those generators supplying power to the electricity system by 
discharging a storage and G'"S vice versa (charging of a storage). The charging 
and discharging efficiency is denoted by ng. Furthermore, all variables might be 
restricted to lower and upper bounds and the starting condition for the storage is 
either a coupling of the first and last time steps or the definition of an initial storage 
level: 


ven < vsr < VO Vs € S,teT 
PI" < per < PR VgeG,teT 
ZO < gi EAN EN teT 


Us,o = Us,NT V Us, = VS VseS 


2.3 Storage Technologies 


Modelling the flexibility provided by actual storage technologies S C S based 
on the equations shown above is straightforward and could easily be applied to 
hydro storages and batteries. Pure hydro pump storages and large-scale battery 
systems allow a rather simple definition of variable bounds. Storage volumes and 
generator capacities are non-negative and time independent, so that the equations 
for the energy level and power output can be reformulated. 


0< vs < VP Ves ter 
0< per < poe Vege GY ET 


O< Pes < PP Vee Gur HET 


While large-scale battery systems allow to ignore external power inflow and 
outflow, which do not result from indirectly connected generators or consumption, 
this is only true for non-cascading pump storage systems. However, we will ignore 
the case of cascading systems in the following and refer to the description of its 
modelling [13]. Including a time series of natural inflow and the restriction of a 
potential minimal outflow, allows for modelling mixed-hydro pump storages and/or 
pure seasonal storages, with an empty set G!”5 in the last case.” 


?In case of missing data for modelling seasonal storages, such as the actual volume or the inflow 
time series, it might be modelled as described in Sect. 2.5. 
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Concerning small scale batteries, such as BESS combined with a PV system, 
the mathematical modelling is identical to large-scale battery systems, as long 
as no further restrictions apply to the battery utilisation. In practical application, 
however, regulators might apply specific restrictions for the operation of such 
systems. In Germany, for example, a widely used support scheme for PV-BESS 
is limiting the grid feed-in up to 50% of the PV capacity.” Given a demand 
Pq, and a maximum grid feed-in BER, these restriction might be easily imple- 
mented: 


D ea — PEM +  pu<) Pas 


eeGPVs gEeGours® deD* 


ys e SPVBESS pe 7 


2.4 Demand Flexibilities 


In general, consumption processes with the ability to store energy can be directly 
modelled as single storages. Under the assumption that a certain share of a specific 
load might be shifted within a certain time range, the storage restrictions, however, 
might be applied to a broader range of demand flexibilities or for modelling the 
flexible potential of aggregated loads. Considering our target to model the impact 
of demand flexibilities at high and extra high voltage levels, we focus on modelling 
the flexibility for large-scale consumption process or of aggregated loads as storages 
SŁ C S. When looking at the load shifting potential of an EV fleet, aggregation 
allows to neglect technical restrictions of single units and to derive the load shifting 
potential from the specific load profile [14]. This finding is supported by [15] where 
an aggregated modelling of BEV as flexible loads led to almost the same results in a 
unit commitment problem as their individual representation, after the BEVs where 
clustered accordingly. 

Therefore, we model the demand flexibility of utility scale HP and the aggre- 
gation of numerous BEV and household heat pumps (HH-HP) by defining certain 
share a of a specific load Py; as flexible pe = af! x P], and adjust the 


general equations based on the load profile in a way such that PÍ e might be 
shifted within a specific time interval. Given a segmentation of the time interval 
T = {TTı,..., TT m} into M subsets, the storage equation might be reduced by 


3KfW-Kredit 275-Erneuerbare Energien-Speicher, from https://www.kfw.de/inlandsfoerderung/ 
Unternehmen/Energie Umwelt/F%C3 %B6rderprodukte/Erneuerbare- Energien- %E2%80%93- 
Speicher-(275) 
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neglecting the variables for storage inflow and discharge: 


t L 
Us t = Us t—1 + > Pest * Ng — Ber VseS”,teT\l 
geGins® 


By restricting the storage outflow to the flexible demand profile 


ZU Se AZO = P ses fer 
deD* 
We define that the flexible demand either directly translates to a system load or is 
stored. In case that load shedding is not allowed, gmin equals Z5“. Assuming that 
the energy of a flexible demand profile has to be consumed within a time interval 
TT m, the bounds of the storage volume are defined by the integral of the profile 
over TT m and equal zeros in the last time step. 


0 if t =max{TT m} 
1 
Ir eTTm deps Pi, else 


Vses!,m=1...M,teT 


max _ _ymin __ 
Vor aa Vst a 


Concerning the bounds for the storage discharge, the minimal or maximal values 
of the profile time series within the certain time interval TT m or the whole time 
horizon can be used. 


. l l 
Pyrex = min{ max a* x pe ‚ max PI!) 


teTT » ? teT 


on I 
pr" = max{ min a * pie , 0} 
i teTTm % 


VdeD’,m=1...M,teT 


In order to model the specific flexible demand technologies, only the definition of 
the time intervals as well as a/ lex, at, a have to be adjusted, where at and a~ 
are scaling parameters within the range of zero and infinity. 


2.5 Renewable Flexibilities 


Given a generation profile, a flexible operation of RES technologies, which depend 
on a storable resource, such as hydro- or bioenergy, might be modelled analogously 
to the case of flexible demand. This might be useful in cases where the actual storage 
parameters are unknown or for modelling virtual power plants, comprising multiple 
small units. The general idea is to model the flexibility by enabling a shift of the 
initial generation profile within a certain time interval. If a certain share aw of the 
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renewable generation is based on stored energy, we first split the generation variable 
into a flexible and inflexible part: 


flex 


min inf lex 
Fe < P + Pst 


max 
8t < Pet 


flex lex max _ pflex 
Por = after x Be: = Por 


Just like in the case of flexible demand, we assume that the flexible generation might 
be shifted within a specific time interval T = {TTı,..., TT m}. In consequence, 
the flexible operation might be expressed by modifying the general storage equation: 


flex in 
Us,t = Us,t—-1 — >: Pg,t /ng + est 


geGouts® 
Ys e SÈ,teT\I 


SR defines the set of storages used to model the generation shift flexibility of 
renewables. By restricting the storage inflow to the flexible generation profile 


min in max __ X flex 
Zur u ost = Zi A Pei 


gEGours® 


We define that the flexible renewable generation is either directly feed-in or stored. 
In case that dumping the flexible renewable generation is not allowed, zmn equals 
Z7". Assuming that the energy of a flexible renewable generation profile has to be 
fed in within a time interval T T m, the bounds of the storage volume are defined by 
the integral of the profile over TT m and equal zero in the last time step. 


ymax _ _ymin _ 0 if t= sodu 
St CA ex 
: Dr ETT m È gegons Pot else 
vse S¥,m=1...M,teT 
Concerning the bounds for the storage discharge, the minimal or maximal values 
of the profile time series within the certain time interval TT m or the whole time 


horizon can be used. 


; l 1 
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; ta l 
P? = max{ min a * Be , 0} 
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In order to model the specific flexible renewable technologies, only the definition of 
the time intervals and a/!@*, at, a7 has to be adjusted. 


2.6 Grid Operation 


Based on the dispatch result for each generation unit and time step from the market 
model, we determine the necessary dispatch adjustment due to grid constraints in 
the second step using a nonlinear optimisation approach. As a conventional cost 
minimising objective function under consideration of grid constraints would lead to 
a nodal pricing result, which would discard the dispatch determined by the market, 
we establish the minimum amount of dispatch adjustment in the entire network 
as the optimisation objective and evaluation figure to assess the efficiency in the 
transmission grid. Therefore, the generalised objective function can be formulated 
as follows: 


. LS 
min ) [Apg] + ) APG 
Ap, Apes 
got dt geG,teT deD 


Positive and negative deviation of generation units from the base solution AP, ; are 
considered equivalently towards the objective, while grid-induced load shedding 
APES is unidirectional by definition. 

We solve the resulting optimisation problem with a multi-period nonlinear 
formulation of the ACOPF utilising an augmented Lagrangian formulation for 
coupling of time-dependent storage units as described in [16]. In order to eliminate 
the possibility of results representing significantly suboptimal solutions caused 
by converging into local optima due to the nonlinear nature of the problem, we 


benchmark our results with a DCOPF formulation, as described in [17]. 


3 Case Study 


In the following, we apply the presented modelling approach on a scenario derived 
from the scenario “Distributed Generation” (DG) of the European Network of 
Transmission System Operators for Electricity (ENTSO-E) Ten-Year Network 
Development Plan (TYNDP) 2016. In our study, we restrict the scope of both 
electricity market and grid to the central European countries and their transmission 
networks, thus not considering the effect of flows from further interconnected 
countries. The considered area contains the countries France, Belgium, Netherlands, 
Switzerland, Austria, Czech Republic, Poland, Germany and Luxembourg. These 
countries form individual market areas, with the exception of Germany and 
Luxembourg, which form a common market area. 


90 M. Ruppert et al. 


In addition to the scenario data from the DG scenario, the closely related high 
renewables scenario C of the german network development plan 2030 was used 
for Germany due to the fact that this source includes data with a higher level of 
detail. Based on the projected national RES capacity development of the scenario 
an optimal allocation planning model as described in [18] was run in the first 
step for regionalisation. In detail, a cost optimal expansion planning with a high 
spatial resolution until 2050 is performed, considering national and regional lower 
and upper bounds for the development of single renewable technologies. Due to 
the longer optimization horizon, exceeding the forecast horizon of the underlying 
scenarios, a minimal RES-E share on the demand is defined in the later years, such 
that an 80% renewable share is achieved in Europe until 2050. 

In the second step, the capacity adequacy is analysed based on a Monte- 
Carlo simulation of the generation availability. Based on an integrated European 
approach, accounting for a cross-border balancing, and the availability of a flexible 
demand shift (BEV and heap pumps in households and district heating) and flexible 
renewables (hydro storage and biomass) we computed the needed dispatchable 
capacity.* Afterwards we utilised results obtained from a generation expansion 
planning problem of the European power system in order to meet the calculated 
capacity demand for a secure operation of the system [19]. 

In the current scenario, we assumed that existing power plants are decommis- 
sioned at the end of their technical lifetime. For coal-fired power plants, a premature 
decommissioning pathway until 2040 is assumed and individually adjusted for each 
unit with regard to its technical parameters and national regulations, aligned as 
much as possible with the political guidelines on coal phase out known today. The 
portfolio of hydro power plants with storage or pump-storage capability is based 
on today’s generation and announced construction or expansion projects under the 
assumption, that generators reaching the end of their lifetime are replaced with the 
same parameters. 

In the presented case study, 555 existing or currently projected and 145 new built 
thermal power plants and 110 hydro power plants are operational in the scenario 
in 2030. The resulting thermal, hydro pump and hydro turbine capacity and the 
location of their connection to the transmission grid are shown in Fig. 1. 

The case study year with input data for transmission grid as well as generation 
and demand is chosen as 2030 and the simulations are performed with an hourly 
time granularity. For the transmission grid part, a weekly optimization horizon of 
transmission grid without rolling horizon was chosen. For this, storage states of each 
unit at start and end time were fixed to the state given by the previously determined 
state after market-based dispatch. 


4The results showed a mismatch between the calculated adequate capacity and the capacity in the 
scenario, which was significantly lower in some countries. The missing DSM modelling might 
explain some part of the mismatch. 
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Fig. 1 Installed capacity per fuel type at the transmission grid substations. The circle volume is 
proportional of the legends circle volume of 2GW 


3.1 Transmission Grid Model 


The transmission grid model contains the interconnected AC network of the same 
countries considered on the market side, including overhead lines and cables of 
voltage levels above 200 kV. For the studied area, this includes the nominal voltage 
levels of 220 and 380 kV. Additionally, high voltage DC (HVDC) lines are included 
in the dataset. HVDC lines connecting offshore wind generation are considered at 
the point of RES-E calculation and allocation and subsequently reduced from the 
dataset. All AC and HVDC lines are connected to busbars, which are assigned to 
georeferenced substations. The number of busbars at each substation is limited to 
the number of voltage levels present at the respective substation. In addition to the 
transformation between the extra high voltage levels, the transformation to the high 
voltage levels between 60 and 150kV is considered and used for the connection of 
small-scale generation and demand as described in [18]. 

The data include the present state of the transmission grid as well as projected 
expansion measures in terms of deconstruction, replacement and construction of 
substations, busbars, lines and transformers until the year 2030. The technical data 
for each grid element was derived from publicly available sources whenever possible 
and otherwise approximated based on available information from predominant 
comparable equipment in the same geographical and organisational area. Among the 
sources used for completing the dataset are the static grid models of the transmission 
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grid operators of the Central Western Europe Region (CWE), as well as the grid 
dataset of the TYNDP 2016 and Open Street Map (OSM). The future expansion 
of the national grid networks was obtained from various sources such as national 
network development plans of the respective grid operators and the projects listed 
in the TYNDP 2016. The model represents a snapshot of the projected grid topology 
in 2030 and includes a total of 3010 busbars, distributed over 1513 substations. The 
busbars are connected by 4103 AC lines and transformers, as well as six HVDC 
lines, with the exception of the interconnector between Belgium and Germany being 
located in the German market area. Furthermore, 237 reactive power compensation 
elements are located in the entire system, with either positive or negative reactive 
power provision. The snapshot of the transmission grid topology for 2030 is shown 
in Fig. 2. 

The interconnection capacities between market areas are determined based on 
the available thermal capacities of interconnecting lines from the transmission grid 
dataset, assuming that 40% of the total thermal capacity of interconnecting AC 
lines between the market areas and the full thermal capacity of interconnection 
HVDC lines is being made available for market operation. Generators from non- 
intermittent sources are assumed to have reactive power provision capabilities of 
cosp = 0.925 w.r.t the installed generator capacity. Furthermore, RES-E with a 
nodal distribution of installed capacity as shown in Fig. 3 is assumed to be able to 


Fig. 2 Map of the transmission 2030 used in the case study. Red denotes 380 kV lines, green 
220 kV lines and purple HVDC lines 
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Fig. 3 Installed renewable capacity per type at the transmission grid substations. The circle 
volume is proportional to the legends circle volume of 2GW 


provide reactive power up to a level of cosp = 0.95 for all sources of renewable 
energy. With regard to the reactive power consumption, the reactive power demand 
in every market area and of every demand type is assumed to remain constant over 
time with a lagging power factor of cosp = 0.95. The maximum voltage deviation 
from nominal voltage is restricted to +10% and grid components can be utilised up 
to their continuous thermal current restriction. 


3.2 Parameters for Flexibility Modelling 


Considering the modelling of new consumers, their flexibility and their regionalisa- 
tion the following assumptions were made: 


3.2.1 Heat Pumps in Households 


The modelling of electricity consumed due to heat demand in households from 
HPs follows the approach of the Munich City Utilities (SMW) for calculating 
the standard load profile of HPs within their grid, based on the yearly household 
load and the temperature. Assuming a full flexibility (a = I, at = 
1, «” = O) within a 3h horizon ([00 : 00 — 03: 00), [03 : 00 — 06: 00),....), 
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the regionalisation follows the same approach as that of the household demand, 
restricted to single and two family buildings. 


3.2.2 Heat Pumps in District Heating 


We apply the same approach for modelling the generation profile and its flexibility 
as it is the case for HH-HP. Due to the availability of a heating grid and the 
possibility to store larger amounts of energy, a full flexibility within a 24h horizon 
[00 : 00 — 24:00) is assumed. The regionalisation is in general based on the 
capacity of combined heat and power (CHP) units within a country. In Germany 
furthermore the data of the German District Heating Association (AGFW) for the 
actual district heating grids demand and location where combined with the CHP 
database for a more detailed analysis. 


3.2.3 BEV PBEV 


Profiles are based on [14, 20] for a direct charging at arrival and adjusted based on 
Table 1. 

The regionalisation corresponds to that of the household demand in all European 
countries except of Germany, where the vehicle registration numbers at NUTS 3 
level were used as a distribution key instead of household number and electricity 
demand. Concerning the flexibility, we assumed the potential to shift the full 
demand (a!!* = 1, at = I, a = 0) within a 12h horizon either during 
the day or the night ([06 : 00 — 18:00), [18 : 00 : —06 : 00). 


3.2.4 PV-BESS 


In the current analysis, decentral BESS installed with PV on household level are 
modelled as simple battery storages and differ only in their sizing and regionaliza- 
tion from classical battery storages. Based on a linear relation between household 
demand and self-optimised PV-BESS capacity, obtained from a mixed integer 
optimization of PV-BESS sizing of households [21] in northern Germany and 


Table 1 BEV parameter set 


Charging at Charging at Mean distance Demand 
home [kW] work [kW] Range [km] | [km/a] [kWh/100 km] 
2020 | 3,5 0 200 13.073 22 
2030 | 5 5 300 13.071 22 
2040 |7 7 400 13.364 22 


2050 |11 11 400 13.481 22 
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southern Germany, an upper bound for the regional distribution of PV-BESS is 
calculated. By adjusting the BESS capacity according to the actual installed rooftop 
PV smaller 15kWp in single and two family buildings, the key for the final PV-BESS 
distribution is computed. 


3.2.5 RES-E Flexibilities 


For the flexibility of hydro storage and biomass generation we assumed a/!¢* 


0.5, at = inf, a = 0, meaning that half of the generation profile is assumed to 
be flexible and restricted to the maximum value of the total time series of the flexible 
part. Furthermore, the potential to shift the generation within 24h was assumed in 
case of biomass and of 168 h in case of hydro storage. 


4 Results 


4.1 Scenarios 


In the following, we will present the results obtained from the simulations performed 
with the input dataset described above. The results are presented for four scenarios, 
which differ due to the amount of utilisation of the flexibilities both on the market 
and the transmission grid side. In this, we divide the available flexibilities into 
traditional flexible generation from hydro storage and into distributed flexibilities 
(HP, BEV, BESS, RES-E flexibilities) In the Base scenario, distributed flexibilities 
remain static according to their inelastic profile and time-dependent generation is 
not available for adjustment during grid operation. Thus, only hydro flexibilities 
are being utilised in the market operation and are subsequently fixed during the 
second step of grid operation simulation. The NoFlex scenario extends the Base 
scenario by adding hydro flexibilities to the available measures while resolving 
transmission grid congestion, keeping distributed flexibilities inelastic. The Flex 
scenario represents a further flexibilisation, with distributed flexibilities included 
in the market operation as described in Sects. 2.2, 2.3, 2.4, and 2.5 Lastly, in the 
DynFlex scenario the available BESS capacities after the market dispatch are added 
to the available adjustments during grid operation, increasing the number of storage 
units considered in the grid operation by 896 nodal distributed elements. 


4.2 The Impact of Distributed Flexibilities 


The yearly aggregated results are shown in Fig.4. Transmission lines which 
reach their continuous thermal current limit during the time horizon are colored 
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Fig. 4 Hours with maximum loading of transmission grid lines and average dispatch adjustment 
per type. The circle volume is proportional to the legends circle volume of 500 MW 


according to the number of hours where this occurs, marking overloaded lines which 
require adjustments in the generation and demand dispatch in order to resolve the 
congestion. Congestions on the transformation level and voltage boundaries are 
omitted from Fig.4. The highest number of fully loaded hours is found on the 
HVDC lines, implying that the start and end locations of these lines are generally 
well suited for relieving the stress on the AC grid. In the AC grid, areas with 
significant congestion can be found in Southern France, the border of France and 
Belgium, Northern Germany and the border area of Poland, the Czech Republic 
and Austria. Generally, the required decrease in generation at certain locations is 
reached by renewable curtailment rather than power decrease of thermal power 
plants. A trend of required power reduction in the southern part of the network 
can be observed over the entire year, while power increase is more focussed in the 
northern and eastern areas. When interpreting the results, possible discrepancies 
between the generation and demand scenario on the one hand and the grid model on 
the other hand have to be considered. For example, the grid expansion considered 
based on the TYNDP in France is significantly lower than the evolution of the grid 
in Germany in the time frame between today and 2030. Thus, the increase of both 
RES-E and new electricity demand types has an higher impact on a system that does 
not undergo a significant transformation. 

The impact of the utilisation of distributed flexibilities on the security of supply 
can be seen in the reduction of the required load shedding in Table 2. While both the 
Base and the NoFlex scenario require load shedding in order to achieve a feasible 
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Table 2 Load shedding results 


Base NoFlex Flex DynFlex 
Hours with LS 240 136 11 9 
Avg. hourly LS [MW] | 20,48 9,00 0,78 0,78 
Max. hourly LS [MW] | 10043 6429,48 1794,14 2156,31 


solution in more than 100 single hours during the year, this value is greatly decreased 
in the Flex scenario and even more in the DynFlex scenario. The main reason for this 
is the reduction in peak demand and generation due to the market-based dispatch of 
flexibilities as described in Sect. 2. Furthermore, the additional control capabilities 
— albeit with a small capacity over a longer time horizon - lead to less congestion 
events, which cannot be prevented without curtailing demand. 

In the Base scenario, the average hourly positive dispatch adjustment requirement 
in the entire transmission grid area considered is 8.672 GW, with peaks reaching 
adjustments of up to 18.5GW. As a permanent additional power provision is 
required in order to balance the transmission losses over lines and transformers, the 
positive adjustment does not decrease to zero even in cases without any violation of 
voltage or thermal limits. In the NoFlex scenario, the hourly requirement reduces by 
48 MW while the peak increases to 20.1 GW. The hourly average decreases further 
to 8.54 GW in the Flex scenario and 8.497 GW in the DynFlex scenario, with the 
peak decreasing to 17.7 and 17.6GW. The average negative dispatch adjustment 
required is not distorted by additional grid losses, over the scenarios it decreases 
by 3.6%, 2.1% and 1.0% from the NoFlex to the DynFlex scenario. Similarly to 
the positive adjustments, the negative yearly peak values are reaching the highest 
absolute values for the NoFlex scenario and the lowest for the DynFlex scenario. 

Overall, the results show that as expected more flexibilities available both in the 
market and transmission grid lead to a lower amount of required redispatch. Yet, the 
different utilisations of both central and distributed flexibilities have a lower impact 
on the grid operation than expected. Among flexibilities the inclusion of large- 
scale hydro in the Flex scenario has the largest impact, as observed in the yearly 
values of negative dispatch adjustment. Operating distributed flexibilities according 
to central market signals leads to a small improvement in terms of transmission grid 
congestion. A similar observation can be made for the inclusion of BESS in the 
adjustable generation on the grid side. The comparably small difference between 
the scenarios might be explained by the discrepancies in the transmission grid model 
and the generation and demand data for 2030, with structural congestion accounting 
for a large part of the observed overloading events. However, the resulting security 
of supply is largely increased by utilising distributed flexibilities with load shedding 
mostly rendered unnecessary, as can be seen in Table 2. In order to improve the 
accuracy of findings for the total amount, more grid expansion in the input data or 
a better coordination of renewable expansion and demand increase on the one hand 
and the transmission grid on the other hand might be needed. 
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4.3 AC vs. DC Formulation 


The comparison of hourly results of the DynFlex scenario with the results of an 
identical simulation performed with the DC restrictions in the load flow equations is 
shown in Fig. 5. For the sake of an equivalent comparison, grid losses are subtracted 
from the positive dispatch adjustment of the AC solution, and the demand of the 
DC model was increased by today’s average national transmission grid losses for all 
market areas considered. Otherwise, both model formulations, restrictions and input 
data were kept consistent. The DC problem was solved using the commercial solver 
Gurobi. Even though the set of DC constraints form a relaxed formulation of the 
problem, the AC solution has a lower requirement for adjustments after correcting 
the AC results for the required adjustments to account for grid losses. This is due 
to the additional positive generation increase needed to account for grid losses in 
the AC case, which are part of the optimisation problem in the formulation chosen 
in this paper and thus contribute to relieving existing congestions. As a result they 
lower the additionally needed adjustment after the subtraction of the losses, leading 
to lower requirements, even though the linearised DC formulation is commonly 
referred to as a relaxation of the AC formulation. In the section up to 5 GW, the 
effect of voltage limits can be seen, which leads to higher AC adjustments as voltage 
is constant in the DC formulation. Overall, individual deviation in single hours 
can be significant while the entire trend and structure are similar. The correlation 
between the positive adjustment value series is 91.1%. This shows that the DC 
approach is suitable as an estimator for the AC equations for the investigation of 
general trends, while individual results and peaks might differ significantly due to 
the simplifications undertaken. For these cases, the AC formulation leads to better 
results. 
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Fig. 5 (a) Distribution of positive adjustments for AC/DC results. (b) Distribution of negative 
adjustments of non-intermittent generators. (c) Distribution of negative adjustments of intermittent 
generators 
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5 Conclusion 


In this paper, we presented a two-step approach to investigate the influence of 
distributed flexibilities on the operation of the transmission grid in the central Euro- 
pean area. For this, we have developed a market optimisation formulation, which 
enables us to obtain dispatch results which include different kind of distributed 
flexible generation and demand sources. Subsequently, we analysed the effect of 
four scenarios with different implementation rates of flexibilities in market and 
grid operation based on a case study using data for the year 2030. Furthermore, 
we compared the results of our multi-period AC formulation to determine grid 
congestion with a linearised DC formulation. 

We conclude based on the required adjustments over the simulation horizon 
of 8760h, that a scenario like the DG scenario of the TYNDP, with high RES- 
E increase and the additional introduction of demand-side flexibilities into the 
electricity system is generally feasible for the anticipated transmission grid of the 
year 2030. Locally concentrated, larger adjustments can be explained by the gap 
between known grid expansion in selected countries and the ambitious targets of 
the chosen scenario. Also, distributed flexibilities show to have a strong effect on 
solving grid congestions that lead to load shedding, which is necessary in the Base 
scenario of our study in more than 200h in order to achieve a feasible solution. This 
is almost entirely resolved when optimising distributed flexibilities with the market 
and enabling the adjustment of BESS during grid operation. 

Furthermore, the results of the four scenarios investigated show that the impact of 
increased participation of both central and distributed flexibilities in the market and 
grid operation has a positive but comparably low influence on the overall required 
adjustment to operate the grid within its technical limits in our chosen case study. 
The largest change on the average adjustments required occurs when adding central 
storage units to the set of generators and consumers available for adjustment in the 
grid simulation. The comparison of the AC and DC results showed that the required 
adjustment can be lower in the AC case, if the provision of grid losses is included 
into the AC formulation and generation increase is performed at locations which are 
suitable for lowering stress on the grid. Additionally, the individual results of each 
hour showed a generally well reproduced trend using the DC relaxation, however 
single hours can differentiate significantly in both cases. 

In the current state a simplified approach of modelling the flexible operation 
of shiftable demands and renewables based on their profile was integrated in the 
market model. In the future we are planning to integrate the decentral flexibilities 
modelling into the grid model und validate it against a more detailed representation 
of individual units. As demand and supply uncertainty have a major impact on the 
utilisation of storages and load shifting potentials, we are furthermore planning 
to extend our modelling to a stochastic optimization. Finally, the extension of the 
current modelling approach to a unit commitment problem would allow us to fully 
evaluate the contribution of flexibilities for the future power system. 
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Part III 
Managing Demand Response 


A Discussion of Mixed Integer Linear A 
Programming Models of Thermostatic Chente; 
Loads in Demand Response 


Carlos Henggeler Antunes, Vahid Rasouli, Maria João Alves, Álvaro Gomes, 
José J. Costa, and Adélio Gaspar 


Abstract In smart grids, it is expected that electricity retailers will offer time- 
differentiated tariffs with significant price variations in short periods. Consumers 
are then encouraged to engage in demand response strategies by making the most 
of their flexibility in the operation of end-use appliances to minimize the electricity 
bill without compromising on comfort. Air conditioning systems are particularly 
suited to be controlled, by profiting from the thermal inertia of building units. This 
paper presents novel mixed integer linear programming formulations to optimize 
the operation of a thermostatically-controlled air conditioning system, thoroughly 
discussing their main features and associated computational difficulties resulting 
from their combinatorial nature. 


Keywords Mixed integer linear programming (MILP) - Thermostatic loads - 
Demand response 


1 Introduction 


Buildings in the commercial and residential sectors account for about 40% of the 
total energy consumption in European Union countries [1, 2]. Heating, Ventilation 
and Air Conditioning (HVAC) systems represent one of the most significant loads 
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contributing to electrical energy consumption. HVAC systems are particularly suited 
to be controlled by making the most of: (i) the thermal inertia of spaces to be 
conditioned (leading to some time dissociation between the energy consumption 
and the energy service provided), and (ii) the flexibility of consumers’ preferences 
regarding the provision of energy services in face of dynamic time-of-use tariffs 
(displaying some willingness to bear the indoor temperature somewhat farther from 
the reference settings for a limited period of time). Since time-differentiated tariff 
schemes are expected to become relevant commercial offers in electricity retail 
markets in smart grids, demand response actions should be developed to achieve 
optimal operation of HVAC systems. Optimal HVAC control can be beneficial for 
consumers (by reducing the energy bill without jeopardizing comfort), retailers 
(by managing buying and selling prices) and grid operators (by contributing 
to release the strain in distribution networks and enhancing the utilization of 
renewable sources). Consumer engagement in demand response programs can be 
made operational by means of automated energy management systems endowed 
with optimization algorithms and the capability to control appliances. Therefore, 
adequate and tractable models for optimal thermostat programming should be 
developed as operational tools (to be embedded with sensors and control) for 
demand response programs. 

The potential of thermostatic loads for demand response actions has been 
exploited in the literature, due to their significant consumption and the possibility 
of being controlled. For this purpose, mathematical models have been proposed, 
in particular mixed integer linear programming (MILP) formulations, i.e. with 
integer and continuous variables. The authors in [3] consider heat pumps for heating 
residential buildings with a floor heating system using a linear state space model in 
an Economic Model Predictive Control framework. Thermal models are developed 
in [4] for a smart house for determining the value of thermal inertia in demand 
response dynamic pricing using a MILP formulation. The authors in [5] present 
a user-centric demand response control for scheduling the electric space heating 
load under price and load uncertainty to minimize a weighted-sum of the expected 
payment, loss of comfort, and financial risk of a customer while considering the 
end-user’s preferences. The household thermal behavior is modeled by means of 
a two-capacity building model. An HVAC system is considered in [6] that is 
controlled (in combination with other type of loads, PV generation and storage) 
by a home energy management system, thus enabling residential consumers to 
participate in demand response programs. A price-based demand response strategy 
for multi-zone office buildings to optimize the energy cost of HVAC units and 
the thermal discomfort of occupants is formulated in [7] as a MILP model. The 
authors in [8] develop an approach based on a partial-differential equation model of 
thermal diffusion to determine the thermostat settings to minimize the electricity 
bill for a consumer with energy time-of-use and power prices, in which the 
optimal thermostat programming for HVAC is formulated as a constrained dynamic 
optimization problem. The authors in [9] consider the optimization of the investment 
and operation planning of a decentralized energy system, which is subject to 
different sources of uncertainties, encompassing photovoltaic generators and load 
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flexibility using heat pumps in combination with thermal energy storage units for 
space heating and domestic hot water, which is tackled by two-stage stochastic 
programming. 

This paper presents novel comprehensive MILP models in which the physical 
characteristics of a thermostatically-controlled air conditioning system are consid- 
ered. These models have the common aim of determining the optimal thermostat 
operation using different forms of thermostat control. It will be discussed that 
the combinatorial nature of these models makes it difficult to solve them by a 
state-of-the-art exact solver. This fact impairs embedding realistic mathematical 
models in automated building energy management systems when a fine-grain time 
discretization of the planning period is considered. 

The study is made in the perspective of the building owner or tenant, i.e. the one 
who aims to minimize the energy bill. These models can also be useful to optimize 
demand response strategies that can be of interest for the grid operator (possibly 
through the mediation of an aggregator of end-users’ flexibility that uses this asset 
to participate in ancillary service provision markets). Thus, this paper adopts an 
Operational Research methodology focusing on the development of mathematical 
programming models. 

The paper begins by developing a simple thermodynamic model to derive a time- 
discretized equation expressing the instantaneous indoor temperature (at instant 
t;) as a function of the indoor temperature, the external temperature and the air 
conditioning system operation at the preceding instant (t;—1). MILP models are 
then presented accounting for the hysteresis behavior of the thermostat between 
minimum and maximum temperatures defining the dead-band around a reference 
temperature (set-point). It is shown that, if just modelling this behavior, the model 
is solved in an acceptable computation time. However, if a minimum ON or 
OFF period is considered (also to avoid excessive switching), offering increased 
operation flexibility, the computation effort becomes impracticable. 

In Sect.2 the building thermal model is developed. MILP models of a 
thermostatically-controlled air conditioning system with different features are 
described in Sect. 3. Section 4 reports extensive computational experiments of the 
different MILP models. Conclusions are drawn in Sect. 5. 


2 Development of the Building Thermal Model 


A space, a building unit, is considered for being heated/cooled by means of an air 
conditioning (AC) system. Assuming the building unit as a (homogeneous) control 
volume, whose instantaneous indoor temperature is uniform and denoted by 0°” (t), 
the overall thermal energy balance on a heat rate basis (at some generic instant fr) is 
[10]: 


doi" 
OE THOE TIMORE RO a (1) 


108 C. H. Antunes et al. 


for AC in heating mode, and 


doi" (t) 


ext g _ AAC zur 
q4” O+ 4 -g M=C FF 


(2) 


for AC in cooling mode. 

In these expressions: g°*'(t) is the instantaneous external load [kW]; q£ (t) is 
the instantaneous internal load (rate of heat generation, by occupants, lighting, 
appliances) [kW]; C = pc,V is the overall thermal capacity [kJ/°C], in which 
p and cp are, respectively, the mass density [kg/m?] and the specific heat at constant 
pressure [kJ/(kg.°C)] of indoor air, if no other thermal inertia of the building is 
considered, or some weighted values for indoors, and V is the volume of the building 
unit [m?]; gee (t) is the instantaneous air conditioning heating/cooling load [kW]. 

Considering a finite-difference approach over a time-step Af to integrate equa- 
tion (1) in a fully explicit time discretization — all variables at the “previous” 
instant t;—1, except for temperature gin (ti) at the “present” instant of Ar - it can 
be assumed that the rate of the energy storage, i.e. the right-hand side in (1) is given 
by cW, and the external load q% (1) © U.A [I (n_1) — 0" (ti-1)]. 
U is the (weighted-average) overall heat transfer coefficient of the building unit 
envelope [kW/(m?.°C)] and A is the surface area of the envelope [m2]. Therefore, 
U.A represents the overall thermal conductance of the unit envelope [kW/°C]. 

Introducing the expressions above in (1), equations (3), (4), (5) are obtained: 


Cosi , 
valent — oi" | + a8, tats = a o, 6) 
. ; U.A ; At 
tar rt int) (4) 
; U.A ; U.A e At 
AAN + a +S) (5) 


In a dynamic simulation of buildings, the internal thermal load q [kW] is usually 
defined according to the operation profile (e.g., daily/weekly profile of occupation, 
of lighting, etc.), i.e. all internal loads (thermal power) associated with the type of 
utilization of the building. 

Without loss of generality for the application of the model, and for the sake of 
simplification in this context, the following assumptions are considered: 


1. the internal heat loads are neglected: a, = 0. 
2.0, SCOP Pes = COP PAE Si 


ici nom 


where C OP is the AC coefficient of performance, Be and pee are the power 


demand and the AC nominal power, respectively, and s;,_, is the AC control variable 
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Table 1 Typical values of parameters a, 6 and y for different discretization intervals of the 
planning period 


Parameter At=1h At = 15min At = 5min At = 1min 
a 0.4276 0.8569 0.9523 0.99046 

B 0.5724 0.1431 0.0477 0.00954 

y PC/kW] 11.1111 2.7778 0.9259 0.1852 


ON/OFF at instant f;_ı. Equation (5) can be rewritten as: 


U.A , U.A At 
g = (1 AN) a + HEAD OF + COP Prom Su (6) 
oi" =a Oi" 4+ BO ty PACs, (7) 
TE (8) 
U.A 

b= (9) 
a) ae (10) 

ars 


Therefore, since COP, PAC, At and C are always positive quantities, the sign 
in equation (10) is positive as this equation has been derived from (1) for heating 
mode. Otherwise, y = -24C O P (negative) for cooling mode. 

The coefficients a, and y derive from the building characteristics (area, 
envelope, etc.) and the CO P of the AC. Considering a building unit with a 225 m? 
floor area and 3 m height, and the properties of the air inside (density of 1.19 kg/m? 
and specific heat of 1.007 kJ/(kg.°C)), the air thermal capacity is C = 810kJ/°C. 
Taking U-values of 0.35 and 0.30 W/(m?.°C) for the building facades and roof, 
respectively, a value of U.A =~ 0.129kW/°C is obtained. Assuming a COP = 
2.5, typical values of parameters a, ß and y are displayed in Table 1 for different 
discretization intervals of the planning period. 


3 Mathematical Models of a Thermostatically-Controlled AC 
System 


This section presents different MILP models aimed at determining the optimal 
operation of the AC within a planning period to minimize energy costs, considering 
distinct forms of control. Although other types of loads could be considered 
(such as shiftable and interruptible appliances), the goal herein is to focus on the 
mathematical modelling of the AC operation. Shiftable loads include dishwashers, 
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laundry machines and clothes driers, i.e. appliances for which the operation cycle 
cannot be interrupted once initiated. The supply of interruptible loads, which include 
electric water heaters and the battery of electric vehicles, can be interrupted within 
preferred time slots provided a certain amount of energy is supplied. In general, 
shiftable and interruptible loads require MILP models, which can be solved very fast 
by state-of-the-art solvers (e.g. Cplex), although requiring many binary variables. 
Examples of such models are the ones presented in [11] and [12]. However, the 
MILP modelling of the thermostat behavior is more complex and imposes a higher 
computational effort to obtain the optimal solution. In particular, MILP models 
are presented accounting for the hysteresis behavior of the thermostat between 
minimum and maximum temperatures defining the dead-band around a reference 
temperature (set-point) established by the user according to his comfort preferences. 

The planning period consists of T time intervals t£ = 1,...,7 of a given 
duration; this discretization can be, for instance, 15-, 5- or I-min. The length of 
the time interval is denoted by Ar, in hours. For instance, if a 1-min discretization 
is used then At = 1/60h. 

Figure 1 displays the behavior of the thermostat of an AC in heating mode. 
The hysteresis operation of the thermostat prevents excessive switching when the 
indoor temperature is around the set-point, which may be established as the mid- 
point within the dead-band defined by the comfort range of indoor temperature 
[orin gmax] specified by the user. These comfort bounds may depend on ¢, i.e. may 
be different during the planning period ([0”"", @/"*],t = 1,..., T); for the sake 
of clarity, we will consider the comfort bounds constant throughout the planning 
period in the formulations presented below. Hysteresis determines the AC control 
variable: 


— sı = 1, i.e. the AC continues ON as the indoor temperature oj” increases above 
gin until it reaches 0”"@* (a in Fig. 1); 

— s, = 0, i.e. the AC continues OFF as the indoor temperature decreases below 
9@* until it reaches 6”""” (b in Fig. 1). 


3.1 Modelling the Thermostat Behavior 


Model 1A forces s; = 1 when the indoor temperature (6! n); 


— is below the lower bound of the comfort band (ein <0") OR 


Fig. 1 Behavior of an AC s=1 oo 
thermostat for heating mode 


S; =0 a 
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— is within the comfort band (0"'" < on < 9") AND in the previous instant 


the AC was also ON (s;_1 = 1) 
Model 1A (heating mode): 


T 
min) [c Ds ‚JAt 
t=1 
s.t.: 
aj” = ao" + BOR + y PAC St- Kerne: 
gi" > gmax _ My, tel... 
ein > grin _ Mw, t= | PENO A 
St-1 < Wr E E 8 
y +w s <1 t=1,...,T7 


s: € {0, 1}, y: € {0, 1}, w € {0,1}, t=1...,T 
a" > 0 V i 


The binary variables s; control the thermostat switching: 


0 if the AC is OFF 
1 ifthe AC is ON 


St 


(1) 


(12) 
(13) 
(14) 
(15) 
(16) 
(17) 


(18) 


The auxiliary binary variables y; and w; establish the consistency conditions 


associated with thermostat switching to the ON position (s; = 1): 


1 if gin < Max 
yt = 


0 otherwise 


1 ifo” <gmin ORs, =1 


0 otherwise 


M is a positive large number, which is used in (13) and (14) to enforce the logical 
conditions. M must be larger than any temperature value; for instance, M = 100 
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(or any higher value) is suitable as temperatures are in °C. The constraints of Model 
1A ensure that: 


1. if oi” < g"in then y; = Land w; = 1 by (13) and (14); hence, by y;+w;—s; < 1 
(16), s; = 1; 

2. if mM < ae < 6" and s;_; = 1, then y; = 1 and w; = 1 by (13) and (15); 
hence, by y; + uw; — st < 1 (16), s; = 1; 

3. in other cases, s; is free to be 0 or 1. The objective function pushes these variables 
to 0, in order to turn the AC to the OFF state and minimize cost. 


Additional input information to be specified include: 


— The temperatures defining the thermostat dead-band, i.e. the minimum temper- 
ature 9”i” (°C) below which the thermostat should be ON (s, = 1) and the 
maximum temperature 0” (°C) above which the thermostat should be OFF 
(s: = 0). 

— The initial indoor temperature 65" (°C) and the initial status (0, 1) of the 
thermostat so. 

— A forecast of external temperatures 0¢*" (°C). 

— The nominal power of the AC, PAC (kW). 

— The cost of energy c; (€/kWh) for t = 1,..., T, in which the coefficients c; 
derive from a time-of-use tariff (a commercial offer of a retailer selected by the 
user). 


This model assumes that the “natural state” of the control variable s; is O (since it 
has a positive coefficient in the objective function to be minimized). The objective 
function forces the binary variables to O whenever the constraints do not impose 
these variables to be 1. However, in this model it may happen that, when indoor 
temperature drops below the upper bound temperature (6’"") after having been 
above it (and hence s, = 0), the control variable switches from O to 1 within 
the thermostat comfort band (note that there are no constraints forcing s, = 0 
because this would be the normal state of the variable). This may be beneficial for 
the objective function by avoiding switching on at some later time intervals when 
the electricity cost is higher. Therefore, Model 1A does not replicate exactly the 
thermostat hysteresis behavior. This model guarantees the maintenance of the ON 
state within the comfort dead-band (a in Fig. 1) but not the maintenance of the OFF 
state (b in Fig. 1). This freedom to switch ON to benefit the cost objective function 
is the reason why this model takes so much computation time. A model forcing the 
control variables to 0 or 1 according to the thermal balance equation and thermostat 
hysteresis, i.e. a rule-based system modelling the thermostat switching, is almost 
instantaneously solved. This is implemented in Model 1B. In addition to Model 
1A, Model 1B (also for heating mode) explicitly forces sı = 0O when the indoor 
temperature CH m); 


— is above the upper bound of the comfort band (6/” > 6’"“*) OR 
— is within the comfort band (6”"” < gi" < gMax) AND in the previous instant 
the AC was also OFF (s;_; = 0) (b in Fig. 1) 
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Model 1B (heating mode): 


min Ct Pros] At (19) 
s.t.: 
o” = a0" + BOM ty PAC sı-ı t=1,... T (20) 
git > grin _ Ms, t=1,...,T (21) 
oi” < om + Mz; t=1,...,T (22) 
gin > gmax _ My, | eee a (23) 
oir < 0” 4 M(l1—s,) t=1,...,T (24) 
zt +y — s1 +s <2 t=1,...,T (25) 
zt + yt +St-1— s < 2 t=1,...,T (26) 
st € {0,1}, ye € {0, 1}, z: € {0,1} +=1,...,T (27) 
ai" > 0 t=1,...,T (28) 


The binary variables s; keep the same meaning as in Model 1A, i.e. they control the 
thermostat switching. 

The auxiliary binary variables y; and z; establish the consistency conditions 
associated with thermostat switching to the ON (s; = 1) and OFF (s; = 0) positions: 


Lae ao 
JS 
0 otherwise 


1 fe? > grin 
Zt = R 
0 otherwise 


Constraints (21) impose that s; = 1 if oj” < 9” Constraints (22), (23) together 
with (25), (26) impose that s; keeps the value of s,_ı within the comfort band. 
Constraints (24) impose that s; = 0 if 6 Pes gmax. 

However, this model is just the mathematical formulation of logical implications, 
which establish the values that the control variable should have according to 
the indoor temperature 0f” and thus no optimization is really at stake. That is, 
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s; variables are rigidly determined by the conditions that establish gi” and the 
thermostat operation (Fig. 1). In this case, there is a single feasible solution (the 
one that complies with the thermostat hysteresis switching rules), so obtaining the 
solution to this model is almost instantaneous. 


3.2 Modelling Minimum ON/OFF Duration Periods 


A new model (Model 2) has been developed to offer the possibility of controlling 
the AC to make the most of time-differentiated tariffs, i.e. being ON during lower 
electricity price periods when it is not strictly necessary to heat the building space so 
that thermal comfort is achieved even limiting the time ON when electricity prices 
are higher. In face of this freedom given to the model, it is necessary to inhibit 
excessive switching, which is accomplished by imposing a new set of constraints. 
In the comfort band, the AC control variables s; are determined to minimize costs 
in face of dynamic tariffs, but once there is a switch from 1 to 0 or from 0 to 
1, then the new state ON/OFF should be kept for at least d time intervals (these 
minimum ON/OFF periods could be different for each state). These new constraints 
establish: 


— (33): if there is a change from sr—1 = 1 to s; = 0, then the subsequent s; should 
be equal to 0 during d time intervals. 

— (34): if there is a change from s;_; = 0 to s; = 1, then the subsequent s; should 
be equal to | during d time intervals. 

- (35) and (36): in the last d time intervals of the planning period, if there is a 
switch from 1 to 0 or from 0 to 1, then this position must be kept until T. 


This model is also quite demanding from the computation time point of view due to 
offering the possibility of switching from | to 0 or from 0 to | in the comfort band. 
Model 2 (heating mode): 


T 
min > [cr pie st] At (29) 
t=1 
s.t.: 
o” = ai", + pom ty PAC si t=1,.. T (30) 
gin > grin _ Ms, Pole (31) 


gi" < gar + M(1 — sr) t= 1,...,T (32) 
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t+d—1 
1-5 +s- +3 2 s; <2 t=1,...,T-d+1 (33) 
J= 
t+d—1 
1-50 +8:-1+7 2 s21 t=1,...,T-d+1 (34) 
T 
les +s a ee t=T—d+2 T (85) 
t t—1 PE RE EEE 
T 
1=s +s iS yai t=T-d+2 T (36) 
t t—1 Dee Te yore, 
s: € {0,1} t=1,...,T (37) 
gi" >0 t=1,...,T (38) 


3.3 Modelling an Inverter AC 


Inverter technology AC appliances have growing acceptance because of increased 
efficiency, extended life and elimination of abrupt load and temperature variations. 
In this type of appliances, a variable-frequency drive adjusts the compressor and 
the cooling/heating output. The previous models can be extended to accommodate 
inverter units, for instance considering they can be operated at different levels with 
respect to the nominal power. This can be modelled by introducing additional binary 
decision variables: 


1 ifthe AC is operating at load level r 
ô = (r =1,..., R) attimet € T 


0 otherwise 


For instance, for R = 5 with the power levels 20-40-60-80-100% of the nominal 
power, the AC operation is modelled in Model 3. In this model, 6”"” and oM«* 
should be interpreted as absolute bounds the user is willing to endure. Note that 
the sets of constraints in the previous models could also have been used here (for 
simplification purposes they were omitted). pe is the power required to operate 
the AC at time ¢ of the planning period, i.e. the AC is either OFF ( pe = 0) or 
supplied at one of the specified power levels. 
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Model 3 (heating mode): 


T 
min 2 [c PAO Ag 
t=1 


s.t.: 


o” =a0 +p ty PAC t=1,... 


pe = (0.25) + 0.457 + 0.68} + 0.8.67 + 57) PAC 


nom 


1 2 3 4 3 

Chae res oe reer ae Seal t=1,... 
git > gmin t=1,.., T 
oi” ge t=1,...,T 

5 € {0,1} r=1,...,5; t=1,... 
pee =0 t=1,...,T 


3.4 Modelling the Cost vs. Comfort Trade-Off 
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(39) 


If the energy price at instant t, c;, is higher than a given threshold, c£, and the user 
is willing to accept the indoor temperature some degrees, 0°, below the minimum 
of the comfort band Hi” (or above the maximum 6” in cooling mode), this can 
be modelled by defining the auxiliary variables y; and the following constraints: 


1 fo > cE 
yt = . 
0 otherwise 


em" << t=1,.. 


cê — ce + My, = 0 f= ye 


c® —c; -M(1— yy) < 0 t=1,... 


BR; (47) 
‚T (48) 
‚T (49) 


The features of the models described above have been presented separately but 
they may be combined. For instance, the different power levels in Model 3 can 
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be combined with the thermostat modelling of Models 1 and 2 and/or with the cost 
vs. comfort trade-off. 

In the next section, illustrative results of the exact solution (whenever possible 
with a given computation budget) of these models using a MIP solver are presented 
drawing attention to the computational effort. 


4 Illustrative Results 


4.1 Data 


A planning period of 24h with a time discretization of | min is adopted, i.e. 1440 
time steps, which imposes a heavier computational burden. The thermostat dead- 
band is defined by 6”""” = 20°C and 0”“* = 24°C. The initial indoor temperature 
a = 18°C in all models, except in model 3 in which it is 65" = 20°C. The 
initial thermostat status so = 0. The external temperatures 9°*'(°C) for a winter 
day (from a meteorological station in Coimbra, central Portugal, on January 1°", 
2012) are available at http://www.uc.pt/en/org/inescc/publications/files/temp_out_ 
winter.xlsx. 

The nominal power of the AC, PAC = 1.5kW. The time-differentiated energy 
prices c;(€/kWh), considering 10 sub-periods, are given in Table 2. Experiments 
have been made on an Intel(R) Core(TM) i5, 3.33 GHz computer with GAMS ver. 
24.0.2 and Cplex ver. 12.5.0.0. 


4.2 Results 


A computational budget of 7200 s or relative MIP gap of 0.5% was established. The 
behavior of thermostat switching as well as the evolution of indoor temperature are 
displayed in Figs. 2 and 3 for Model 1A and Model 1B, respectively. The cost of 
the corresponding optimal solutions for Model 1A and Model 1B, as well as the 
dimension of the models and the Cplex time to obtain these solutions are given 
in Table 3. A minimum duration ON/OFF after a thermostat switching occurs is 
considered with d = 3 and d = 5min (Model 2). The behavior of thermostat 
switching and the evolution of indoor temperature are displayed in Figs. 4 and 5 for 
Model 2 with d = 3 and d = 5 min, respectively. The cost of the corresponding 
optimal solutions ford = 3min and d = Smin, as well as the dimension of the 
model and the Cplex time to obtain these solutions are given in Table 4. Considering 
that the AC can operate at power levels 20-40-60-80-100% of PAC in Model 3, 
the cost of the optimal solution as well as the dimension of the model and the Cplex 
time to obtain these solutions are given in Table 5. The corresponding solution of 
the AC operation and indoor temperature variation is displayed in Fig. 6. 
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Fig. 2 Optimal AC operation and indoor temperature for Model 1A 
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Fig. 3 Optimal AC operation and indoor temperature for Model 1B 
Table 3 Optimal solutions, problem dimension and Cplex time (Model 1A and 1B) 
Optimal Relative |Cplex |N. binary | N. continuous 
Model name solution (€) | gap (%) | time (s) | variables | variables N. constraints 


Model 1A 2.416125 |43.47 |7200 14320 1440 7199 
Model 1B 2.465125 |0 10.01 14320 1440 10078 
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Fig. 4 Optimal AC operation and indoor temperature for Model 2 with d = 3 min 


The frequency of switching in Model 1A is higher than in Model 1B due to 
the freedom given by that model to minimize cost. In Model 1A, the AC is ON 
for a longer time in pricing sub-period P3 to account for the increase in price 
in subsequent sub-periods (Fig.2). Note that the solution presented in Fig.2 and 
Table 3 for Model 1A still presents a MIP gap of 43.47% after 2h of computation, 
due to the combinatorial difficulties of this model. Model 1B, as expected, is solved 
to optimality instantly. 

The temperature peaks in Figs. 4, 5 and 6 are due to the capability of the models 
to accommodate for comfort by making the most of lower price sub-periods. The 
MIP gaps in Model 2 (Table 4) are still significant after 2h of computation. 

The optimal solution to Model 3 is better than the solutions obtained for the other 
models due to the possibility of the inverter AC to work at different power levels. 
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Fig. 5 Optimal AC operation and indoor temperature for Model 2 with d = 5 min 
Table 4 Optimal solutions, problem dimension and Cplex time (Model 2) 
Optimal Relative} Cplex | N. binary} N. continuous 
Model name solution (€) | gap (%) | time (s)| variables | variables N. constraints 


Model 2 (d = 3 min) | 2.015875 7200 | 1440 14320 8636 
Model 2 (d = 5 min) | 2.07725 10.62 |7200 | 1440 [4320 8636 


Table 5 Optimal solutions, problem dimension and Cplex time (Model 3) 


Optimal Relative Cplex N. binary |N. continuous 
solution (€) | gap (%) time (s) | variables variables N. constraints 


2880 7200 


Power(kW) 


a SM mind ala ot manmanm 


Time (minute) 


Fig. 6 Optimal AC operation and indoor temperature for Model 3 


5 Conclusion 


This paper presented a set of mixed integer linear programming models in order 
to optimize the operation of a thermostatically-controlled air conditioning system 
(AC), aiming at minimizing the energy costs. These models capture different 
physical control characteristics: 


— Model 1A models the hysteresis operation of the thermostat (in heating mode), 
preventing from excessive switching when the indoor temperature is around a set- 
point. This model gives, however, some freedom to thermostat switching when 
the indoor temperature is within the dead-band, enabling to change from OFF to 
ON if this benefits the cost objective function. 

— Model IB removes the flexibility given by the previous model, forcing the AC to 
keep its state when the indoor temperature is within the dead-band. This model 
implements, in a MILP formulation, a rule-based system of logical implications. 
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The feasible region is reduced to one solution and the model can be solved very 
quickly, unlike Model 1A which is computationally very demanding. 

— Model 2 gives some flexibility while inhibiting excessive switching, by consid- 
ering that, when the indoor temperature is within the dead-band, the thermostat 
must keep its state (ON or OFF) for at least a predefined number of time intervals. 
This model is also quite computationally demanding. 

— Model 3 introduces a new feature that can be used to extend the previous models: 
it models inverter units, considering that the AC can operate at different levels of 
nominal power. 


Modelling cost reduction at the expense of accepting worsening the indoor temper- 
ature by some degrees was also introduced in this paper. It has been shown that the 
combinatorial nature of these models imposes a severe computation burden, which 
in some cases impairs obtaining the optimal solutions in a practical, acceptable 
computational time by a state-of-the-art solver. Therefore, this work lays the 
foundations for understanding those computational difficulties and gives insights for 
the development of other approaches, namely dedicated meta-heuristics customized 
for the features of the models, in order to offer good solutions with a suitable 
computational effort. 
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Weighted Fair Queuing as a Scheduling A 


Algorithm for Deferrable Loads in Smart wes 
Grids 


Tuncer Haslak 


Abstract Weighted Fair Queuing (WFQ) is implemented for the problem of load 
scheduling in demand side management. Power demand, wait time and group-to- 
group fairness are the basis for three variants of the algorithm’s implementation. The 
results are compared to a Greedy strategy with regard to the residual of renewable 
power supply and the suggested measures of fairness. WFQ proves comparable 
to Greedy in terms of the primary objective and in addition is capable of equally 
distributing resources and distress caused by deferral. 


Keywords Demand side management - Optimization - Weighted fair queuing 


1 Introduction 


Renewable energies are a resource that strain the grid through intermittent availabil- 
ity and difficulty in prediction. Demand side management addresses the issue by 
reversing the paradigm of grid operation and controlling power consumers instead 
of only power generators. 

Within demand side management there is a need for robust algorithms that face 
the unpredictability aspect of renewable energies, which makes real-time algorithms 
with no forecast information favorable. 

Finding solutions that improve the residual power problem is only the first 
step. Said problem is the question of how unused renewable supply should be 
handled. Using the Greedy algorithm significant improvements can be accom- 
plished. However, the results indicate that distress would be unevenly spread in the 
consumer population. If distress is distributed unevenly, compensation would also 
be distributed unevenly, which means that the market is not design towards fairness. 
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This paper introduces Weighted Fair Queuing (WFQ) as a scheduling algorithm 
for deferrable loads in smart grids. Three variant implementations are presented 
that feature different concepts of fairness: serving more and simultaneously smaller 
loads, fairly distributing wait times, and treating any number of groups equally. 


2 Materials and Methods 


2.1 Model 


The experimental environment is designed in AnyLogic and uses an interface to 
control power consuming processes, as presented in [1]. The interface is suitable 
for any process that is capable of prolonging periods of its activity or inactivity as 
can be seen in Fig. 1. This requires no intimate knowledge of the process: when 
an algorithmic decision is made to alter the behavior of an individual load, the 
request is passed on to the agent representing the load. In accordance with its 
internal conditions it then accepts or rejects. This is to replicate the sovereignty 
of loads, as in actuality algorithms are informed of the current states and consider 
them accordingly. 

Every process is divided into 4 states: activity and inactivity, and for each of 
those one deferrable and one non-deferrable. Processes can alter their power demand 
from step to step or follow a mathematical function—anything that describes the real 
behavior. The model captures distinct periods in which switching can be deferred. 


Fig. 1 State automaton that governs the state changes for deferrable loads; States can be changed 
by progression of time (clock) or request (question mark) 
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Fig. 2 This graph shows the power demand of a generic electric vehicle over approximately 


2.5h with a peak power demand of 22kW. The simulation contains 300 consumers, this graph 
exemplifies one specimen 
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Fig. 3 The automated welding machine requires 150kW over 1-3 s, repeated 12 times over. In 
the final step once 230 kW are needed for 5s. This graph serves to show the internal working of a 
consumer agent 


The deferral and operation times are a matter of survey in the companies (interviews, 
data sheets and surveillance). The total load shape is based on real loads that are 
replicated with the simulation. Figures 2 and 3 are two examples of the power 
demand and the repeated activation of loads. 
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The analyzed algorithms make no use of forecasting of any kind; they operate 
based on information of the “now”. The power supply is processed data-point-by- 
data-point, left to right. 


2.2 Population 


The model population consists of approximately 300 individual loads. These encom- 
pass beverage production, metal working, glass finishing, plastics production, textile 
treatment, sewage treatment, electric vehicles and compressors for both pressure 
applications and freezing. The switching behavior of a load for the purposes of the 
simulation is detailed. The realistic accuracy was verified by a test implementation 
in [2]. The composition of the population, i.e. the proportions of types of loads 
are elected to reflect a midsized city by approximation, based on statistics found 
in [3, 4] and [5]. The overall peak power demand across the whole population is 
25 MW, with single members requiring power between 1 and 350 kW. The average 
population power demand is 5.8 kW (cf. Fig. 7). State changes are deferrable by 15% 
to 50% and depend on the current and future state (arrows in Fig. 1) of the individual 
load. At any given time an average of approximately 20% of loads is deferrable. 


2.3 Objective Function 


Considering the negative impact of renewable energies, optimizing the load popula- 
tion toward a renewable energy supply profile seems adequate. As the volatility of 
the renewables is the straining property, immediate consumption by the deferrable 
loads would lead to the minimization of the residual power function and thereby 
stabilization. Residual power (R) is the amount of power leftover from renewable 
supply after subtracting the demand. A sample from historic supply data is picked, 
that consists of wind and solar energy over 24h with a resolution of one data point 
per 3s. This deviates from the resolution of simulated demand, which is continuous. 
Figure 4 shows the sample. At this resolution the objective function incorporates no 
useful gradient information. 

However, selecting a random supply function is not sufficient: First the data 
points of the supply function must be scaled according to the demand of the popula- 
tion of deferrable loads. The first simulation experiment is a simple observation of 
how much power the consumer population requires when there is no optimization at 
all. This is the so called “natural run” that allows any and all requests to switch on or 
off. This unimpeded 24 h simulation returns a power demand function (P (t) Demand) 
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Fig. 4 Sample supply function, sum of wind and solar energy in northern Germany 


in Watts. The integral thereof is the energy demand (Wpemand) in Watt-hours 
(Eq. 1). 


24h 
/ P(t) Demand dt = WDemand (1) 
0 


In Eq. (2) the integrals of the natural run and the supply functions which are 
WDemand and Wsupply are equated. This yields a result for c, the scaling factor. The 
result of this approach is a function of power supply (P(t) supply) that can satisfy 
the demand in terms of energy (kWh), but in terms of the progression of power over 
time (kW) any algorithm must manage transient deficits. 


! 
WDemand = €: Wsupply 
P(t) Demand # Pt) Supply 


(2) 


This scenario seems manageable concerning not only algorithms, but is also 
realistic with regard to resource allocation. The latter meaning that purely renewable 
supply seems unrealistic and would probably feature a safety factor to account 
for fluctuations in demand. At the same time, however, in an economic sense a 
rough balance between supply and demand is sensible, because of investment and 
operative costs. At the least, equal supply and demand appears to be a good initial 
metric. 

To facilitate the comparison the summation from equation (3) is used, approxi- 
mating the integral of the residual power for discrete time steps: 


T 


R= > | Psupply — Paemana|(t) (3) 
t=0 


This condenses the result of one 24h simulation into one manageable value. 
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2.4 Algorithms 
2.4.1 Greedy 


The Greedy strategy (Algorithm 1) is a heuristic search algorithm. It is incomplete, 
but fast and very simple. As [6] shows it generally yields good results, though it 
is not capable of systematically finding global optima. I chose this algorithm as a 
starting point for its popularity, non-specificity and low computational requirements. 
I have previously outlined the use of the Greedy algorithm in [1] and [2]. The 
algorithm sorts the available candidates according to size and attempts to activate as 
many as possible, as can be seen in Algorithm 1. 


2.4.2 Weighted Fair Queuing 


Weighted Fair Queuing (WFQ) is a sophisticated algorithm designed by Demers 
et al. (cf. [7]), building on the work of Nagle [8]. WFQ manages congestion in 
gateway queuing of datagram networks. The problem occurs when participants 
in a network attempt to send data packets, but some entities always send larger 
packets than others. WFQ replaces the pragmatic FIFO (first in, first out) paradigm 
by categorically distinguishing datagram sources by packet size and permanently 
assigning them to queues with predetermined weights. Weight, packet size and a 
queues history determine which packet gets serviced next. 

Identifiable parallels to the scheduling of deferrable loads lead to the application 
of WFQ in a novel way. The deferrable loads also feature a well-defined requirement 
for resources which is their power demand in Watts. This parameter is equatable 
to the packet size, as is the network bandwidth to the usable residual power. The 
significant reason that makes a more sophisticated algorithm desirable is the unfair 
dissemination of resource with the Greedy strategy. By sorting according to size, 
Greedy categorically prefers large loads, which disadvantages smaller loads. 


Algorithm 1: Greedy 


1 while forever do 


2 forall Load do 
3 if Residual > Load(i) then 
4 activate; 
5 calculate Residual; 
6 else 
7 | deactivate; 
8 end 
9 end 
10 end 
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WFO (2) employs a heuristic called virtual finish time! F (Eq.4), which is 
calculated using the packet size, or in this case the wattage of a load / and the 
weight w of each queue 7. Every load is usually permanently assigned to a queue, 
as are weights to queues. 


F; =— (4) 


This yields a numeric value for each load that wants to schedule for activation. 
For each queue a so-called session virtual finish time SF; is calculated (Eq. 5). 


i—1 
SRh=F+) Fy (5) 
n=0 


Within every queue, all elements are queued on a FIFO basis. The element of the 
queue i to be scheduled next is F;. All previously successfully activated elements 
are added with the summation term F,,. At every point in time the algorithm selects 
the queue with the lowest SF; value. The significant difference in comparison to the 
original WFQ implementation is the explicit check (Algorithm 2, In. 10) whether 
the current residual power is sufficient for the regarded load. In the network realm 
decreasing bandwidth will not cause a packet to be rejected service-it will, however, 
take increasing amounts of times to submit. This “flow” property is not given in 
the case of deferrable loads, therefore “choking” cannot be applied, and switching 
decisions are discrete. 


Fairness 


In datagram networks fairness is the property of allotting a certain resource to any 
participant according to a predefined key. Such a key maybe simple and assign 
equal fractions to everyone. In the case of WFQ, however the distribution scheme 
accounts for a priori observations in terms of the typical resource requirements of 
each participant. With WFQ bandwidth cannot be manipulated which means that the 
amount of sheer data transmitted will at best remain the same (though more packets 
cause more gaps). In fact the objective is to service a greater number of participants. 


Fairness of Wattage 


This exact thought is transferred to demand side management applied to a load’s 
wattage. In contrast to the original network application participants do not forfeit 


‘heuristic, not physical time 
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Algorithm 2: Weighted fair queuing 


1 begin function calculateSFs 
forall queues do 

Fi) <Il/w; 

SF(i) < F(i) + F(n); 


end 
while —Abort do 


2 
3 
4 
5 end 
6 
7 
8 calculate Residual; 


9 calculateSFs forall Load do 
10 if Residual > Load(i) \ Restart then 
11 activate; 
12 Restart < true; 
13 else 
14 | Abort < true; 
15 end 
16 end 
17 end 


after multiple trials in vain. It is implausible to assume that a load be denied access 
altogether. Loads forcing themselves to power on is unfavorable, but in line with 
the conventional grid operation paradigm of supplementation in such a case. As 
the key influenceable property is the length of part of the “on” and “off” episodes 
of each load, disproportionately large loads compensate with their “on” time. This 
means that larger loads are operated for shorter periods of time by the algorithm. 
The objective of employing this queuing scheme is to serve more participants in the 
same time interval. Here fairness is targeting decongestion. 


Fairness of Distress 


“Wait time” is a second, additional fairness measure worthy of investigation. 
Wattage is a strictly objective parameter. If we consider the period that a process 
can wait to be served, this requires very intimate knowledge of the process to 
determine the validity or honesty of time parameters as declared by a machine 
operator. A company might easily misrepresent this information. For the simulation 
experiment at hand the processes were personally surveyed and documented during 
live observations. Especially as there was no motivation presented to falsify data, 
the time parameters can considered to be sufficiently accurate. 

In this regard fairness would disregard the wattage of a load and consider the 
time it has been waiting to get served. Fairness in this sense is the equal distribution 
of waiting-distress among the entirety of the load population. 

For this variant fixed queue assignments are replaced with dynamic ones. The 
longer a process waits, the further it moves into queues that make it more likely to 
be picked next. To this end standardized wait time is favorable, i.e. the fraction of 
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time waited in proportion to the maximum amount of time the load can wait (Eq. 6). 
Without this measure loads that provide less flexibility are irrationally rewarded. 


twaited 
Is = —— (6) 


tmax 


Fairness Group-by-Group 


As for a third measure of fairness the load population is subdivided into regions of 
equal size, representative of quarters or districts. In this configuration each queue 
represents one region, and regions compete for the available resources. This option 
is designed to grant each group the same right to resources. Any attribute can be used 
to queue in this manner. This generic queuing scheme targets an equal distribution 
of deferral times and power to independent groups, but not within the groups. As 
the virtual finish time paradigm is still applied, smaller loads are preferred. 


Weights and Queues 


There is no definitive answer on how weights should be chosen. Existing procedures 
such as that outlined in [9] are not applicable because the prevailing problem of node 
congestion in networks is not a relevant phenomenon. Therefore a mode of selecting 
weights should rely on the available parameters. Each of the fairness measures 
suggested above requires its own weighing scheme. 


2.4.3 Weight by Power (WFQ-Power) 


Considering that the resource “power” is engrossed correlating with wattage, a 
load’s wattage should inversely determine its likelihood of being elected. As a 
number of queues to which all loads are assigned, 5 was elected as it is typical. 
The fifth of the population with the highest wattage is assigned to the according 
quintile, and so forth. The weight for each queue is determined by computing its 
average wattage. Following equation (4), fairness is established when the virtual 
finish times of all individual loads are equal. This means solving for w and equating 
the average wattages of all queues. Table 1 shows which weights are assigned to 
queues, and the respective average wattage within the queue. 

For example: At 250kW load # 1 is a member of the queue weighted at 0.09, 
i.e. the group of highest wattages. Compared to bottom quintile member load # 2 of 
20kW and its weight of 1. Load # | has a virtual finish time 138 times higher. The 
mismatch in wattage reflects proportionately in the weights. 
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Table 1 Averages of power 


: , u of power in queue | Queue weight 
in queues and corresponding 


weights (WFQ-Power) 220 kW 0.09 
120kW 0.17 
90 kW 0.23 
40kW | 0.5 
20kW 1 


Table 2 Loads for WFQ-Time are dynamically assigned to queues by the percentage of wait time 
expended from their maximum so far. In order to augment urgency, queue weights are doubled 
from queue to queue 


Expended wait time Queue weight 
80-100% 1 
60-80% 0.5 
40-60% 0.25 
20-40% 0.125 
0-20% 0.0625 


2.4.4 Weight by Wait Time (WFQ-Time) 


Queuing by standardized wait time is carried out dynamically. In contrast to 
wattage-queuing there is no a priori assignment. With increasing wait time, a load is 
assigned to queues with higher weights that equate to higher likelihoods of getting 
picked. Once again there are five queues, each designated for a 20% interval of 
wait time completed. Table 2 shows how weights are governed for each queue. 
Decreasing weights augment urgency. 


2.4.5 By Region (WFQ-Geo) 


In the case of grouping by region equal weights are assigned to all queues. As 
the fictitious regions comprise equal portions of every category of the overall 
population, and given the equal weights, this queuing scheme acts comparable to 
a round-robin. This is true for all types of WFO with equal weights (cf. [8]). 


3 Model Limitations 


The proposed model requires highly automated processes. The underlying work 
plans are machine accessible so that the ensuing steps and their duration can 
be projected. Semi-automation or increased human interaction would hinder this 
procedure. 

Even though they were not encountered during survey, processes that can 
continuously adapt the level of their power input can only be replicated as a series of 
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discrete steps. On the contrary this model is a response to the lack of continuously 
adaptable processes. Yet this is a limitation. 

In addition the size of the consumer population cannot be arbitrarily increased. 
To this end an arbitration mechanism appears feasible. Multiple optimization algo- 
rithms would work on a fraction of a larger problem which in turn is consolidated 
by optimization layers higher up in the hierarchy. This approach could also help to 
reduce the amount of datagrams that are sent between the loads and the optimizer. 


4 Results 


In each simulation there is a population of approximately 300 manageable loads that 
must be scheduled over 24 h following a randomly selected supply signal. There are 
five experiments: 


1. Natural run: All loads operate without any algorithmic interference. This is a 
reference trial. 

. Greedy 

. WFQ: Power 

. WFQ: Wait time 

. WFQ: Region 


nb Wh 


Results are compared regarding the fulfillment of the primary optimization 
objective from equation (3)-the minimization of residual power deviations. The 
three queuing schemes of the WFQ variants have different secondary objectives. 
The algorithm in experiment 3, with weights defined by power demand, aims to 
serve more participants, especially smaller ones. Experiment 4 entails queuing by 
wait time, which is designed to equally distribute the wait time percentages over the 
entire population. Queuing by region, experiment 5, is designed to allocate the same 
resources to arbitrarily defined subsets of the population. 

In order to judge the fulfillment of these secondary objectives additional param- 
eters must be analyzed. Cycles counts all activations which is the sum of every 
individual power-on that was granted-in the datagram sense this is analogous to 
serviced packages. The standardized average wait time is indicative of how long 
every activation was deferred in percent of the maximum. Deferral length is the 
sum of all deferrals as time. The variance of this value is representative of how 
evenly deferrals are distributed among the population. 

Table 3 summarizes the results of experiments 1-4. As outlined in Sect. 2.3 the 
significant data is the power consumption function that every algorithm causes. 
Figure 5 is the consumption profile caused by the Greedy strategy. Consumption 
curves for all algorithms including the natural run (experiment 1) can be found in the 
appendix (Figs. 7, 8, 9, 10, and 11). The visual differences in the residual curves are 
exiguous due to the algorithmic similarity. The key difference is the residual energy. 
Systematic differences do not express discernibly in the graphs. Showing which load 
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Table 3 While Greedy performs best in terms of residual, it causes a high variance in deferral 
lengths (inequal waiting). WFQ-Power and WFQ-Time accomplish their secondary objectives of 
equal distress distribution (variance) and increasing the serviced loads (cycles) at the cost of energy 


optimization (residual) 

WFQ-Geo 
Residual [MWh] 6749.9 6720.7 
Avg. standardized wait [%] |O 0 21340 1347 | 34.2 35.9 
Avg. deferral length [h] o [383 ET 3.90 3.93 
Variance of deferral length jo [02027 [0.1730 | 0.0973 0.169 
Cycles (switch on & off) 439,278 | 443,366 


8 000+ : : : : 


Power [kWh] 


0 ; + 
00:00 04:00 08:00 12:00 16:00 20:00 24:00 
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Fig. 5 Power consumption for the Greedy algorithm, residual R in upper right corner 


Table 4 This table shows the distribution of wait times of all queues of the WFQ-Geo algorithm. 
The objective here is to treat all queues equally, which hinders global optimization objectives. 
Fairness reflects in the similarity of waiting times 


Queue E op e r 
Avg. standardized wait [%] 35.8 36.0 
Avg. deferral length [h] 3.89 3.92 
Variance of deferral length 0.168 0.164 
Cycles 88,450 |87,921 


is activated at any given time offers no insights from algorithm to algorithm, as no 
pattern emerges. 

The Greedy algorithm performs best by reducing the residual by 54.8%. WFQ- 
Time improves the deferral distribution by a factor of 2.08, while maintaining a very 
similar deferral length. 

Table 4 shows the results of all queues of experiment 5, each representing 
one region. As this algorithm is specified not to improve the global population, 
but to serve queues equally, its results cannot be directly compared to the other 
experiments. The only exception is the residual. Queue deferral length is on average 
11% higher compared to WFQ-Power. Queue by queue comparison shows that 
results group close together for all parameters, but worse in general. 


Weighted Fair Queuing for Deferrable Loads 


135 


Greedy 54,8% 
WFQ-Power 
WFQ-Time 51,6% 
WFQ-Geo 51,1% 
49,0% 51,0% 52,0% 53,0% 54,0% 55,0% 56,0% 
© Improvement on Residual 
Greedy 69,3% 
WFQ-Power 86,2% 
WFQ-Time 88,3% 
WFQ-Geo 84,0% 
0,0% 20,0% 40,0% 60,0% 80,0% 100,0% 
E 80% of distress distributed on [%] of population 
Greedy 
WFQ-Power 84,7% 
WFQ-Time 
WFQ-Geo 76,7% 
65,0% 70,0% 75,0% 80,0% 85,0% 90,0% 


E Load demand satisfaction 


Fig. 6 The WFQ algorithms can compete with the versatile Greedy algorithm in terms of the 
residual optimization (blue). However Greedy focuses 80% of distress on only 69% of the 
population (green). In addition WFQ-Power rejects especially few loads (red) 


Figure 6 summarizes the advantages of WFQ: it is an improvement on the Greedy 
algorithm which is fast and versatile. At the cost of slightly worse results with regard 
to the optimization of the power residual, WFQ algorithms distribute distress more 
equally, and serve more loads all else being equal. 


5 Discussion 


5.1 Improvement 


As to be expected the natural run that features no algorithm or adaptation to 
the objective function causes a high residual. The Greedy algorithm delivers the 
greatest improvement in terms of residual. At the same time this algorithm unevenly 
distributes the wait time distress across the population which is expressed in the 
variance of the average deferral time. WFQ-Time significantly narrows the variance 
of deferral, following its secondary objective. The WF Q-Power variant is capable of 
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placing more individual activations which can be seen in its cycle count. Keeping in 
mind that the objective function cannot be fully attained, as indicated in Sect. 2.3, 
the algorithm still reduces the gap in serviced packages by approximately 60% in 
comparison to Greedy. 

The results from the algorithm WFO-Geo in experiment 5 are inferior in general 
as it causes longer wait times and the variance thereof indicates that wait times are 
heterogeneously distributed within queues. The objective of equal queue treatment 
is achieved at the cost of worse results all in all. 


5.2 Distress 


Distress distribution is a desirable trait in demand side management that is not 
trivially deducible. Wait time is a parameter that cannot be easily subverted, if 
for example the maximum deferral time is repaid as an incentive. Weighted Fair 
Queuing is capable of delivering comparable result to a Greedy strategy in the 
primary objective of residual power improvement. In addition it can be adapted 
to pursue additional objectives. I presented three possible measures of fairness: 
equality between multiple queues or regions (WFQ-Geo), the even distribution of 
wait times (WFO-Time), and serving more participants (WFO-Power)-the latter 
being most similar to the original algorithm. This shows that WFQ can be used 
to introduce fairness to the selection of deferrable loads. 


5.3 Methodological Contribution & Policy Advice 


Despite its complexity Weighted Fair Queuing is a staple of network congestion 
management. It is highly functioning and ubiquitous because of its performance and 
its ability to accomplish fairness. The application to energy distribution problems is 
novel. At the same time it introduces the problem of fairness into the matter, which 
was disregarded because of the focus on residual power optimization, monetary 
compensation or welfare gain. 

The optimization problem presented and solved in this paper is significant in 
its size regarding multiple dimensions. Firstly the time resolution exceeds the 
commonly quoted 15 min intervals that stem from the tertiary grid balancing realm. 
Grid stabilization however can only be achieved when the shrinking number of 
spinning masses can be replaced. The realization that compensatory measure must 
move below the 30-seconds-threshold is at the core of this research. Secondly, 
the population is substantial with 300 members and 25 MW peak. Each load is 
replicated as an agent with internal processes, variables and decision making. 
Furthermore their states are not estimated or stochastically approximated—decision 
making and transmission is acute, meaning that statuses are inquired, evaluated, 
decided on and requests dispatched. Thirdly, no simplifying assumptions are applied 
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to the target function (power supply) such as smoothing. The complexity of the 
problem is embraced to the extent that the resolution of the objective function 
(supply) is the restricting factor. The argument is in favor of preparedness, as 
methodologically the question answered here is that demand side management 
features the necessary scalability and agility. 

The suggested algorithms are implementations of centralized algorithms. As 
the technology is still at an early stage centralization is acceptable. Especially as 
the methodological exploration of this topic-which this research is part of-is still 
ongoing. For future applications, however, this is not advisable, as a single entity 
will possess all information on all clients (loads). The controlling entity would be 
vulnerable to attacks . The presented model, is acting on a request basis which means 
that a load can reject any suggestion to change its state. Designing intrinsically safe 
systems increases complexity, but is quintessential policy advice. 


5.4 Practical Implications 


For companies owning deferrable loads the practical implications are low, especially 
if processes are highly automated. The test implementation presented in [1] required 
no manual interaction. While loads were in interference the optimization paused and 
resumed automatically. The more a deferrable load is interwoven into the process 
of a company, the more difficulty can arise from automated unexpected deferrals. 
Deferrable loads are an effort that aim to improve the usage of renewable loads. 
Some resulting distress is to be expected. Although not insinuated here, market 
design should compensate for this. 

The suggested interface is at the core of the research presented, as it allows for 
the load to negotiate deferral and reject it. From a load’s perspective WFO-Time is 
the most favorable, as it guarantees that no load is overburdened with deferral. The 
Greedy algorithm is unfavorable, as some loads are asked to defer significantly more 
often than others (unfairness). 


5.5 Algorithm Comparison 


There is a multitude of optimization algorithms that can be used. The advantage of 
the above presented algorithms is their simplicity. They can be categorized as local 
search algorithms, which means that they use a type of heuristic, like SF; in Eq. (5). 
As there are no simplifying assumptions, and every load is designed to be addressed 
with a direct request, the problem presents itself to be massive and without useful 
gradient information, which excludes classic optimization algorithms. 

A suitable alternative would for example be “Simulated Annealing” (cf. [10]). 
It is an algorithm that begins by amply exploring the solution space with random 
solutions. A solution in this case is an allowable sequence of activation and pause 
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for all loads. The algorithm continues by rejecting or accepting results, based on 
the metaheuristic property (“temperature”). In principle this algorithm converges on 
the global optimum solution, which local search algorithms are not systematically 
capable of. The significant disadvantage is considerably higher run time. A prelimi- 
nary test implementation shows that, after 14h Simulated Annealing delivers results, 
which are comparable to the Greedy strategy which takes approximately 15 min on 
the same desktop computer (2.5 GHz Intel Core i7 with 16GB 1600 MHz DDR3 
memory). 


6 Related Work 


Gellings laid out the fundamentals of demand side management in [11]. 

The objective function and its definition are based on [12] and [13]. 

Weighted Fair Queuing and the approach of Generalized Processor Sharing are 
based on [7, 8] and [14]. These original sources were used for the adaptation to 
demand side management. 

Most related work is based on [15] where an optimal deferrable load control 
problem is defined. The focus is placed on market parameters which is why price 
bounds are employed. [16] focuses on the problem laid out in this publication. To 
solve this problem a decentralized algorithm is employed, which communicates the 
power residual to all participants. This is in contrast to centralized solution efforts, 
but capable of scheduling multiple instances of the same load. [17] utilizes a similar 
approach but only for singular placement of loads in 24h. 

[18] uses a stochastic model which defers one activation of each load and no 
reactivation. The underlying idea is very strongly connected to the paradigm of 
grid operation. This means that the time scales of balancing power in the European 
grid are adopted and all activity takes places in time steps of 15 min, and deferral 
times are at least 4h. The significant difference is the lack of understanding and 
anticipating load behavior. The publication focuses on the simulation of wind power 
only and the cost of demand side management. 

[19] has a similar approach in that only one activation instance of every 
participant is considered. The multitude of loads is approached as a combinatorial 
problem. Here the smallest time step is one hour. In neither of these loads can extend 
their operative time. 

Formulated as a constraint problem regarding the availability of the demand 
side, [20] manages a day-ahead scenario by interacting with energy markets. By the 
introduction of welfare they move away from monetary quantification. They outline 
the difference in time scale between markets (larger steps) and the load behavior 
(smaller steps). 
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Appendix 


Fig. 7 Power consumption for a non-optimized reference run; 1h equals 100 simulation units 
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Fig. 9 Power consumption for the WFQ-Power algorithm 
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Fig. 10 Power consumption for the WFQ-Time algorithm 
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Fig. 11 Power consumption for the WFQ-Geo algorithm 
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Planning and Operation of Distribution 
Grids 


Cost Optimal Design of Zero Emission A 
Neighborhoods’ (ZENs) Energy System eed 


Model Presentation and Case Study on Evenstad 


Dimitri Pinel, Magnus Korpas, and Karen B. Lindberg 


Abstract Zero Emission Neighborhoods (ZEN) is a concept studied in particular 
in the research center on ZEN in smart cities in Norway to reduce the CO» 
emissions of neighborhoods. One question coming along this concept is how to 
design the energy system of such neighborhoods to fit the ZEN definition[1]. From 
this definition we extract the C O2 balance, requiring an annual net zero emission 
of CO: in the lifetime of the neighborhood. This paper proposes a MILP model 
for obtaining cost optimal design of ZEN’s energy system and demonstrates it on a 
case study. Different technologies are included as investment options and, notably 
PV as a mean of producing electricity on-site. Wind turbines are not included in this 
study because they would not be suitable in the context of most cities. The results 
highlight the importance of PV investment in reaching the ZEN requirements. For 
example, around 850 kW of solar is needed for our test cases of 10,000 m? of floor 
area, for an annual energy demand of around 700 MWh of electricity and 620 MWh 
of heat. The investments in other technologies are small in comparison. 


Keywords ZEN - Sustainable neighborhoods - Zero emission Neighborhoods - 
Energy system - C O2 emissions - Optimization 


1 Introduction 


A ZEN is a neighborhood that has a net zero emission of C O2 over its lifetime. 
Many aspects are embedded in the idea of ZEN. Energy efficiency, materials, users 
behavior, energy system integration are all aspects that need to be accounted for 


D. Pinel (4) - M. Korpås 

Deparment of Electric Power Engineering, Norwegian University of Science and Technology, 
Trondheim, Norway 

e-mail: dimitri.q.a.pinel@ntnu.no; magnus.korpas @ ntnu.no 


K. B. Lindberg 
SINTEF Community, Oslo, Norway 
e-mail: karen.lindberg @sintef.no 


© The Author(s) 2020 145 
V. Bertsch et al. (eds.), Advances in Energy System Optimization, 
Trends in Mathematics, https://doi.org/10.1007/978-3-030-32157-4_9 


146 D. Pinel et al. 


in this concept. In addition, different parts of the life cycle can be included but in 
this paper we only consider the operation phase and no embedded emissions. Two 
types of action exist to make neighborhoods more sustainable. One is to act on 
the demand, via better insulation, user behavior or other efficiency measures. The 
other is to act on the supply and have a local energy system minimizing the C O2 
emissions. There is consequently a need for a way of designing the energy system 
of such neighborhoods. The questions to be answered are, which technologies 
are needed to satisfy the demand of heat and electricity of a neighborhood, and 
how much of it should be installed so that it is as inexpensive as possible. The 
problem is then to minimize the cost of investment and operation in the energy 
system of a neighborhood so that it fulfills the ZEN criteria. This paper presents 
an optimization model to solve such problems with a focus on operations research 
methodology. 


2 State of the Art and Contribution 


The ZEN concept is specific to this particular project, however similar topics have 
been studied in different settings either at the neighborhood level, the city level 
or the building level, for example during the research center on Zero Emission 
Building. In this context, K B Lindberg studied the investment in Zero Carbon 
Buildings [2] and Zero Energy Buildings [3] which are variations around the concept 
of Zero Emission Buildings. In both papers an optimization based approach is used 
to study the impact of different constraints on the resulting design. The second 
one [3] in particular uses binary variables to have a more realistic representation 
of the operation part (part load limitation and import/export). In [4], Gabrielli et 
al. tackle the problem of investment and operation of a neighborhood system and 
show an approach allowing to model the system complexity while keeping a low 
number of binary variables. It also constrains the total CO» emissions. It uses 
design days and proposes two methods for allowing to model seasonal storages 
while keeping the model complexity and reducing the run time. In [5], Hawkes 
and Leach look at the design and unit commitment of generators and storage in a 
microgrid context using 12 representative days per season in a linear program. It is 
particular in that it defines how much the microgrid would be required to operate 
islanded from the main grid and include this in the optimization. It also discusses 
the problematic of market models within microgrids. In [6], Weber and Shah present 
a mixed integer linear programming tool to invest and operate a district with a 
focus on cost, carbon emission and resilience of supply. A specificity of this tool 
is that it also designs the layout of the heat distribution network taking into account 
the needs of the buildings and the layout of each area. It uses the example of a 
town in the United Kingdom for its case study. In [7], Mehleri et al. study the 
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optimal design of distributed energy generation in the case of small neighborhoods 
and test the proposed solution on a Greek case. Emphasis is put on the different 
layouts of the decentralized heating network. In [8], Schwarz et al. present a model 
to optimize the investment and the energy system of a residential quarter, using 
a two stage stochastic MILP. It emphasizes on how it tackles the stochasticity of 
the problem in the different stages, from raw data to the input of the optimization, 
and on the computational performance and scalability of the proposed method. In 
[9], he also studies the impact of different grid tariffs on the design of the system 
and on the self-consumption of the PV production. In [10], Li et al. separate the 
investment and the operation into a master and a follower problem. The master 
problem uses a genetic algorithm to find the optimal investment while a MILP is 
used to find the operation in the follower problem. In [11], Wang et al. also use a 
genetic algorithm, but at the building level and using a multi objective approach 
focused on environmental considerations. A life cycle analysis methodology as 
well as exergy consumption are used to assess the design alternatives. In [12], 
Mashayekh et al. uses a MILP for sizing and placement of distributed generation 
using a MILP approach including linearized AC-power flow equations. In [13], 
Yang et al. also use a MILP approach for the placement and sizing problem but 
consider discrete investment in technologies at the district scale. These papers give 
us an overview of different methods for optimal investment in the energy system of 
neighborhoods or buildings, but none apply the ZEN concept and the influence of 
tight requirements on the C O2 emissions on the modelling and on the results has 
not been demonstrated. 

In this paper, the focus is put on getting a fast yet precise solution that can take 
long term trends, such as cost reduction of technologies or climate. To this end, the 
proposed model uses a full year representation, ensuring a correct representation 
of seasonal storage of heat and electricity, and allows to divide the lifetime of the 
neighborhood into several periods, each represented by one year. It is also different 
by using the Zero Emission framework on a neighborhood level as a guide for 
the emission reduction constraint. This adds an integral constraints coupling each 
timestep and increasing the complexity of the problem. The use of binary variables 
is limited to the minimum. 


3 ZENIT Model Description 


ZENIT stands for Zero Emission Neighborhoods Investment Tool. It is a linear 
optimization program written in Python and using Gurobi as a solver. It minimizes 
the cost of investing and operating the energy system of a ZEN using periods, with 
a representative year in each period. Different technologies are available, both for 
heat and for electricity. It is most suited for greenfield investment planning but can 
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also take into account an existing energy system. The objective function is presented 
below: 
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The objective is to minimize the cost of investing in the energy system as well as 
its operation cost. 

The operation phase can be separated in different periods during the lifetime of 
the neighborhood, and one year with hourly time-steps is used for each period. In 
addition to technologies producing heat or electricity, there is also the possibility 
to invest in a heating grid represented by the binary bhg that also gives access to 
another set of technologies that would be inappropriate at the building level. In the 
equation above, the € represent discount factors either global for the whole study (3) 
or for each period (2). They are calculated in the following way: 
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The calculation assumes that reinvestment in this technology is made for the 
whole lifetime of the neighborhood, and is discounted to year 0. The salvage value 
is also accounted for. The formula used is: 
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In the objective function, Yp P represent the total export from the neighborhood. 
It is simply the sum of all exporty from the neighborhood: Yt, p 
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The most important constraint, and what makes the specificity of the “Zero 
Emission” concept, is the C O2 balance constraint. It is a net zero emission constraint 
of CO2 over a year. It takes into account the emissions from the used fuels and 
electricity with the corresponding C O2 factors for the emission part and the exports 
of electricity for the compensation part. In this study the same factor is used for 
imports and for exports of electricity. This constraint is expressed below, Y p: 
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In the particular ZEN framework of this study, the idea behind the compensation 
is that the electricity exported to the national grid from on-site renewable sources 
allows to reduce the national production, and thus to prevent some emissions from 
happening. The corresponding savings, the compensation, stand on the right-hand 
side of the equation. In the ZEN framework, this constraint is set as an annual 
constraint. It can however also be used for shorter periods of time. 

Other necessary constraints are the different electricity and heat balances which 
guarantee that the different loads are served at all times. The electricity balance 
is represented graphically in Fig. 1. The corresponding equations are also written 
below. The electricity balance is particular because, we want to keep track of the 
origin of the electricity sent to the battery. It is managed by representing each battery 
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Fig. 1 Graphical representation of the electricity balance in the optimization 
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as a combination of two other batteries: one is linked to the on-site production 
technologies, while the other is connected to the grid. It allows to keep track of 
the self-consumption and to differentiate between the origin of the energy for the 
CO; balance. 

Node I (8) represents the main electric balance equation while II (9) and V (10) 
are only related to the on-site production of electricity. Node II (9) describes that the 
electricity produced on-site is either sold to the grid, used directly or stored, while 
node V (10) states that at a given time step what is stored in the batteries is equal to 
what is in excess from the on-site production. 

Electricity balance I: Yt, p 


imp gbi dch 5 „pb-selfe selfc 
Yr,p p oP Doi pest 7 Yt, p,est ee Nest + > 8g,t,p 
8 


est 


(8) 
= dent 2) dinsne + 2, Enep: 4 
b hp 
Electricity balance II: Vt, p, g 
l 
&g.t,p = =; Tr gele T en g (9) 


Electricity balance V: Yt, p 


h b_ch 
eg ype (10) 
& 


est 


Heat also has its own balance, that guarantees that the demand of each building 
is met: 


5 drap +9,9 arno 


yEeQ\HP hp 


dch ch 
+ ) Nhst ` dt, p,hst = X Hb,t,p - Ab + dt,p 
hst b 


(1) 


Note that the demand is not divided between domestic hot water (DHW) and space 
heating (SH). 

The batteries are represented, as mentioned earlier and as seen on Fig. 1, as two 
entities: one on the on-site production side and the other on the grid side. This 
means that we have two “virtual” batteries with their own set of constraints as well 
as constraints linking the two. 
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The first constraint is a “reservoir” type of constraint and it represents the energy 
stored in the battery at each time-step: Yt € T*, p, est 


pb __ „pb pb_ch pb_exp pb_selfc 

Ur pest = Vi—1,p,est + Nest * Yi-1,p.est ~~ Yi-1,p,est zu Yi-1,p,est (12) 
gb _ gb gb_imp gb_exp gb_dch 

Ur pest U], p,est + Nest ` Yi-1,p,est ~~ Yi-1,p,est = Yt—1,p,est (13) 


Equations (14), (16) and (17) link both batteries. They make sure the sum of the 
stored energy in the “virtual” batteries is less than the installed capacity, and making 
sure the rate of charge and discharge of the battery is not violated. Yt, p, est 


pb gb bat 
Vi pest + Ur. pest S Yt, pest (14) 
bat < (15) 

Ur p,est = Xbat,est 

pb_ch gb_imp bat 
Yt, pest $ Yi, pest Z Ynax,est (16) 

gb_dch gb_exp bat 
Yt, pest + Yr,p.est Z Ynax.est (17) 


The storage level at the beginning and the end of the periods should be equal. 
Vp, est 
b b 
ken = ee pesi (18) 


The heat storage technologies also have the same kind of equations as the 
batteries, for example: Yt € 7%, p, hst 


heatstor __ „‚heatstor heatstor 


ch dch 
Ur p,hst = U _],p,hst + N st "Ir,p,hst — Ut, p,hst (19) 


Equations (14) to (18) also have equivalents for the heat storages. However the heat 
storages are not separated in two virtual entities since there is no export of heat from 
the building. 

The power exchanges with the grid are limited depending on the size of the 
connection: Yt, p 


(ye? + yi? FY aa T < GC (20) 


est 
In order to not add additional variables, the mutual exclusivity of import and export 


is not explicitly stated. It is still met however due to the price difference associated 
with importing and exporting electricity. 
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In addition to the above equations, different constraints are used to represent the 
different technologies included. The maximum investment possible is limited for 
each technology. Vi: 

ee (21) 


The amount of heat or electricity produced is also limited by the installed capacity: 


Yq, t, P:4q.t.p Š Xq (22) Yg, t, P : 8g,t,p Š Xg (23) 


The amount of fuel used depends on the amount of energy provided and on the 
efficiency of the technology: respectively Yy € F N Q, p, t and Yy € EN Q, p,t 


dy.t, Ay,t, 
ae (24) dyip = Z (25) 
Ny Ny 


For CHPs technologies, the Heat to Power ratio is used to set the production of 
electricity based on the production of heat. Vt, p 


GdCHP\,t,p 
8CHP,t,p = — (26) 
ACHP 


For the heat pumps, the electricity consumption is based on the coefficient of 
performance (COP). 
Vhp, b,t, p 


Ghp,b,t,p 
dnp.b,t,p = COPhp.b.t.p (27) 

The heat pumps are treated differently from the other technologies because they 
are not aggregated for the whole neighborhood but are separated for each building. 
This is because the COP depends on the temperature to supply, which is different 
in passive buildings and in older buildings and which is also different for DHW and 
for SH, and dependent on the temperature of the source. The source is either the 
ground or the ambient air depending on the type of heat pump. The COP is then 
calculated using a second order polynomial regression of manufacturers data [3] 
and the temperature of the source and of the outside timeseries. The possibility to 
invest in insulation to reduce the demand and improve the COP of heat pumps is not 
considered. The global COP is calculated as the weighted average of the COP for 
DHW and SH. 
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The solar technologies, solar thermal and PV, also have their own set of specific 
constraints. Vt, p: 


Bip tert =, -xpv IRR: p (28) 
IRR 
ST t,p 
go) = xsr a 2 29) 
ap Gstc 


The hourly efficiency of the PV system is calculated based on [14], and accounts 
for the outside temperature and the irradiance. This irradiance on a tilted surface is 
derived from the irradiance on a horizontal plane that is most often available from 
measurements sites by using the geometrical properties of the system: azimuth and 
elevation of the sun and tilt angle and orientation of the panels. 

The irradiance on the horizontal plane data comes from ground measurements 
from a station close to the studied neighborhood which can for example be obtained 
from Agrometeorology Norway.! The elevation and azimuth of the sun is retrieved 
from an online tool.” This calculation takes into account the tilt of the solar panel 
and its orientation. Several assumptions were necessary to use this formula. Indeed, 
the solar irradiance is made up of a direct and a diffuse part and only the direct part 
of the irradiance is affected by the tilt and orientation. However there is no good 
source of irradiance data that provides a distinct measurement for direct and diffuse 
parts in Norway as far as the authors know. Thus we make assumptions that allow us 
to use the complete irradiance in the formula. We assume that most of the irradiance 
is direct during the day and that most is diffuse when the sun is below a certain 
elevation or certain azimuths. This assumption gives a good representation of the 
morning irradiances while still accounting for the tilt and orientation of the panel 
during the day. On the other hand, this representation overestimates the irradiance 
during cloudy days, when it is mostly indirect irradiance. Obtaining direct and 
diffuse irradiance data would solve this problem. 


4 Implementation 


The model presented in the previous section has been implemented in the case of 
campus Evenstad, which is a pilot project in the ZEN research center [15]. This 
implementation of the model and the parameters used are presented in this section. 
Campus Evenstad is a university college located in southern Norway and is made 
up of around 12 buildings for a total of about 10,000 m?. Most of the buildings 
were built between 1960 and 1990 but others stand out. In particular two small 
buildings were built in the nineteenth century and the campus also features two 


!Imt.nibio.no 
Sun Earth Tools: https://www.sunearthtools.com/dp/tools/pos_sun.php 
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Table 1 Technologies used Inv. Cost | Life- Efficiency 
in the Evenstad case and their ; 
main parameters (EKW) | time (%) 
Technology (Years) 
Building | 
PV 1600 25 18 
Solar Thermal 0 |25 70 
Air source HP 556 15 COP, 
Ground source HP | 444 15 COP, 
Biomass Boiler 350 20 85 
Electric Boiler 750 30 100 
Gas Boiler 120 25 95 
Neighborhood 
Gas CHP 739 25 45th; 35e1 
Biomass CHP 3300 25 40ih; 25er 
Heat Pump 660 25 COP, 
Electric Boiler 150 20 100 
Gas Boiler 60 25 [95 


recent buildings with passive standards. The campus was already a pilot project in 
the previous ZEB center and one of those buildings was built as a Zero Emission 
Building. In addition, on the heating side a 100kW CHP plant (40kW electric) 
and a 350kW Bio Boiler both using wood chips were installed along with 100 m? 
of solar collectors, 10,000 L of storage tank, 11,600L of buffer tank and a heating 
grid. On the electric side, the same CHP is contributing to the on-site generation as 
well as a 60 kW photovoltaic system. A battery system is already planned to be built 
accounting for between 200 and 300 kWh. Based on this we assume in the study an 
existing capacity of 250 kWh. We keep those technology in the energy system of the 
neighborhood for one part of the study. In addition, the heating grid is kept in all 
cases (bng = 1). 

The technologies included in the study are listed in Table 1 along with the 
appropriate parameters. 

Two main sources for the parameters and cost of the technologies are used as 
references for the study. Most of the technologies’ data is based on a report made 
by the Danish TSO energinet and the Danish Energy Agency [16] on technology 
data for energy plants. The other source includes the technology data sheets made 
by IEA ETSAP [17] and is used in particular for the gas and the biomass CHP. The 
cost of PV is based on a report from IRENA [18]. The two efficiencies reported 
for the CHP plants correspond to the thermal and electrical efficiency, noted by a 
subscript (;n for thermal and .; for electrical). Note that: at the neighborhood level, 
only ground source heat pump is considered (Table 2) and that PV is only considered 
at the building level but the roof area limit to the size of the PV is not implemented. 

The heat storage values are based on a data sheet by ETSAP [17] while the values 
used for the batteries are based on a report from IRENA [19]. 
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Table 2 Storage 
technologies used in the 
Evenstad case and their main 


Inv. Cost | Lifetime | Efficiency 
Technology | (€/kWh) | (Years) |(%) 


parameters Battery 350 15 94 
Heat Storage | 75 20 | 95 
Table 3 Fuel cost and C Op Fuel Cost (€/kWh) | CO; Factor (gC O>/kWh) 
an Gas 0.055 277 
Biomass 0.041 WE; 
Electricity pee 17 


The values in Table 3 come from different sources. The cost of biomass comes 
from EA Energy Analyses [20], the cost of gas is based on the cost of gas for 
non household consumers in Sweden? (we assume similar costs in Norway). For 
the technologies in Table 1, the O&M costs, expressed as a percentage of the 
investment costs, are respectively: 1, 1.3, 1, 1.3, 2, 0.8, 2.3, 4, 5.5, 1, 1 and 5. For 
the storage technologies in Table 2, the operating cost is 0. The C O2 factors of gas 
and electricity for Norway are based on a report from Adapt Consulting [21] and 
the C O2 factor for biomass is based on [22]. 

The electricity prices for Norway are based on the hourly spot prices for the Oslo 
region in 2017 from Nordpool.* On top of the spot prices, a small retailer fee and the 
grid charges are added.” The prices are rather constant with a fair amount of peaks 
in the winter and some dips in the summer. This cost structure is close to the actual 
structure of the electricity price seen by consumers. We assume hourly billing due 
to its relevance to prosumers and its emergence in Norway. 

The irradiance on the horizontal plane and temperatures are obtained and used 
in the calculations as described in the previous section. The ground station used to 
retrieve data is Favang, situated 50 km to the west of Evenstad. The electric and heat 
load profiles for the campus are derived from [23]. The load profiles are based on 
the result of the statistical approach used in these papers and the ground floor area 
of each type of building on the campus. In addition, the domestic hot water (DHW) 
and Space Heating (SH) are derived from the heat load based on profiles from a 
passive building in Finland where both are known [24]. 

The problems are solved on a Windows 10 laptop with a dual-core CPU (i7- 
7600U) at 2.8 GHz and 16 GB of RAM. Each case typically has about 450,000 rows, 
600,000 columns and 2,400,000 non-zeros. They are solved using the barrier method 
in Gurobi in about 150s each. 


3http://ec.europa.eu/eurostat/statistics-explained/index.php?title=File:Gas_prices_for_non- 
household_consumers,_second_half_2017_(EUR_per_kWh).png 


*https://www.nordpoolgroup.com/Market-datal/Dayahead/Area-Prices/ALL1/Hourly/?view= 
chart 
Shttps://www.nve.no/energy-market-and-regulation/network-regulation/network-tariffs/statistics- 
on-distribution-network-tariffs/ 
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5 Results 


The optimization was run several times with different conditions. It was run with 
a yearly CO» balance with and without including the energy system that already 
exists at Evenstad. When the pre-existing energy system is included, the pre- 
existing amounts of heat storage, PV, solar thermal and biomass heating (CHP and 
boilers) represent the minimum possible investments in those technologies for the 
optimization. The energy systems resulting from those optimizations are presented 
on Fig. 2. 

Both cases are interesting. Indeed the case with the pre-existing technologies 
included in the optimization allows to know in which technology to invest to move 
towards being a ZEN for the campus Evenstad while the case that does not include 
the pre-existing technologies allows to see how it would look like if it was built today 
from the ground up using the optimization model presented here and the given ZEN 
restrictions. 

A first observation from Fig.2 is that the technologies already installed (heat 
storage ST, biomass boiler BB, CHP, battery) do not get additional investments, 
except for PV which gets a lot of additional investments to meet the ZEN criteria. In 
addition to the large investment in PV the only additional investment for Evenstad 
appears to be a heat pump. In the case without any pre-installed technologies the 
system is quite different. There is still a need for investment in PV, though it is 
slightly lower and the optimization does not chose to invest in a battery. On the 
heating part the chosen design uses heat pumps and electric boiler in addition to a 
heat storage smaller than already installed in Evenstad. 
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Fig. 2 Resulting energy system 
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Fig. 3 Self consumed electricity (blue) and total consumption (red) of electricity in the ZEN. (a) 
Summer. (b) Winter 
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Fig. 4 Import (red) and Export (blue) of Electricity from the ZEN. (a) Summer. (b) Winter 


The results highlight the predominance of PV in the results. This shows that the 
other possible designs are not cost competitive. Alternative designs, for example 
relying on biomass CHP, could be incentivized to obtain a better mix of technology. 
The amount of the incentive could be explored by a sensitivity analysis using this 
model, but this remains as future work. 

On Fig. 3, the self consumption and the total demand of electricity is presented 
while on Fig.4 it is the imports (red) and exports (blue) of electricity that are 
presented. Both figures show a week for the case of the yearly balance and including 
pre-existing technologies. In the summer the neighborhood produces electricity in 
excess and needs to send it to the grid. The battery, that is part of the pre-existing 
technologies, is used but is not large enough to allow for relying on self produced 
electricity during the night. It is also not large enough to limit the amount of 
electricity sent to the grid. Figure 4a illustrates this: the exports during the days 
have highs peaks that represent around four times the night imports in terms of 
peak power. This has implications on the sizing of the connection to the grid and is 
especially important in the context of the introduction of new tariffs based on peak 
power in Norway. Indeed, the introduction of smart-meters enables the use of more 
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complex grid tariff structures. Such tariffs would promote avoiding large peaks in 
consumption. This may be beneficial to highly flexible neighborhoods such as ZENs 
and might promote investment in batteries. Investigating the impact of grid tariffs on 
the design of ZENs remains as future work. Outside of the ZEN context, a positive 
impact of certain grid tariff designs has been shown on self-consumption and peak 
electricity import [9]. In the winter, some of the electricity is still self consumed due 
to the CHP that is part of the pre-existing technologies. This self consumption stays 
limited and no electricity is exported. 

Ultimately, all resulting designs require huge investment in PV to attain the status 
of ZEN. In those systems, which rely heavily on electricity, heat pumps and electric 
boilers appear to be the preferred heating solution. 


6 Limitations 


This study has several limitations, on the methodology and on the case study. For 
the case study, assumptions were necessary due to the lack of data, in particular 
for the loads or the insolation (diffuse and direct). For the methodology, the will 
to limit the use of binary variables meant leaving out constraints such as part 
load limitations which would be needed to have a better representation of some 
technologies. In addition, using an hourly resolution leads to an underestimation 
of the storages and possibly of the heating technologies size. There is a trade off 
between the solving time and the precision of the results and the resolution needs 
to be chosen accordingly. Additionally, being deterministic, the model leaves out 
several uncertainties. Those uncertainties concern the evolution of the price of 
the technologies, the electricity price or the price of other fuels and the climate 
conditions. Those can be partially addressed by specifying additional periods in the 
model. The short-term uncertainties are not included either and induce an overly 
optimistic operation of the system. Despite those limitations it provides insights in 
the design methodology that can be used to design the energy system of a ZEN. The 
choice of C O3 factors for electricity is also greatly impacting the results and this 
should be studied in more detail in future work. 


7 Conclusion 


This paper presented in detail the ZENIT model for investment in Zero Emission 
Neighborhoods as well as its implementation and the results on a realistic case study 
of campus Evenstad in Norway, with a focus on methods from the field of operations 
research. The model is formulated as a MILP, using as few binaries as possible. 
The Zero Emission constraint complexifies the problematic of designing the energy 
system of a neighborhood and the long term trends can be accounted for by defining 
periods. For Evenstad, the results suggest that additional investments, mainly in 
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PV, are necessary in order to attain the status of ZEN. Investments happen at both 
levels but mainly at the building level. When the technologies already installed at 
Evenstad are not included, they are not invested in (except for heat storage). The 
optimal choice in order to become Zero Emission for Evenstad in the current ZEN 
framework thus appears to be a massive investment in PV and a heating system 
fueled by electricity. Further work includes disaggregating the heat part of the model 
and a more detailed operation part in the optimization. 

There are key takeaways for policy makers in this study, in particular for 
Norway due to the setting of the case study. The results suggest that the Zero 
Emission constraint used in this study is sufficient to get PV investment without any 
additional incentive. However, under the C O2 factor assumptions used in this study, 
huge investment in PV are made which would be problematic in case of a large- 
scale application of the concept of ZEN. This suggests the need for incentives in 
alternative technologies such as biomass CHPs in case the concept of ZEN becomes 
more common. The methodology presented in this paper can be used to assess such 
policies and their potential effect on investments in ZEN. Other policies such as the 
grid tariff structure can also be studied with this model. Finally, the hourly limitation 
on electricity export from prosumers has recently been replaced in Norway by a 
tariff on exported electricity. The results of this paper suggest that without this 
change in policy, ZEN would become even more expensive due to the necessity of 
large batteries to make the exports more constant. We thus recommend continuing 
on the path of facilitating the development of the number of prosumers for example 
with the implementation of capacity grid tariffs. For countries other than Norway, 
similar methodology can be used to assess the cost and design of ZEN. Further 
policy recommendations cannot be drawn from this study due to the specificity of 
the Norwegian electricity mix, that is reflected in its electricity C O2 factor. 
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Nomenclature 


Indexes (Sets) 


t(T) Timestep in hour within year € [0, 8759] 

b(B) Building type 

yr Year within period € [1, YR] 

p Period 

i(Z) Energy technologies, J = FUEUHPUSU OST VEST; T = 
OUG: 
fa) Technology consuming fuel (gas, biomass, . . .) 
e(E) Technology consuming electricity 


hp(HP) Heat pumps technologies 
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q(Q) 
8(9) 


Parameters 


disc 
Ci 


Variables 
Xi 


ffp 
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s(S) Solar technologies € ST, PV 
qst(QST) Heat storage technologies 
est(EST) Electricity storage technologies 


Technologies producing heat 
Technologies producing electricity 


Discounted investment cost, technology i with re-investments and 
salvage value [€/kWh] 

Discount factor, period p with discount rate r 

Duration of the study [yr]: D = P x YR 

Number of periods in the study [-] 

Annual maintenance cost [% of inv. cost] 

Price of fuel of technology g, period p [€/kWh] 

Electricity spot price [€/kWh] 

Electricity grid tariff, period p [€/kWh] 

Retailer tariff on electricity, period p [€/kWh] 

Charge/Discharge efficiency of battery est [-] 

C O3 factor of electricity [g/kWh] 

C O2 factor of fuel f [g/kWh] 

Heat to power ratio of the CHP [-] 

Size of the neighborhood grid connection [kW] 

Maximum possible installed capacity of technology i [kW] 
Electric specific load of building b in timestep t in period p 
[kWh/m?] 

Aggregated area of building b in the neighborhood [m?] 

Heat specific load of building b in timestep t in period p 
[kWh/m?] 

Efficiency of technology i [-] 

Coefficient of performance of heat pump hp in building b in 
timestep t in period p [-] 

Maximum charge/dis- rate of battery [kWh/h] 

Maximum charge/discharge rate of heat storage [kWh/h] 
Efficiency of the solar panel in timestep t in period p [-] 

Lifetime of technology i [yr] 

Cost associated with a heating grid for the neighborhood [€] 


Capacity of technology i: fori € {f U e U h U s} [kW]; fori € 
{qst Uest} [kWh] 
Fuel consumed by technology f in hour t [kWh] 
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Abstract Grid-connected battery storage systems on megawatt-scale play an 
important role for the integration of renewable energies into electricity markets 
and grids. In reality, these systems consist of several batteries and inverters, which 
have a lower energy conversion efficiency in partial load operation. In renewable 
energy sources (RES) applications, however, battery systems are often operated at 
low power. The modularity of grid-connected battery storage systems thus allows 
improving system efficiency during operation. This contribution aims at quantifying 
the effect of segmenting the system into multiple battery-inverter subsystems on 
reducing operating losses. The analysis is based on a mixed-integer linear program 
that determines the system operation by minimizing operating losses. The analysis 
shows that systems with high modularity can meet a given schedule with lower 
losses. Increasing modularity from one to 32 subsystems can reduce operating 
losses by almost 40%. As the number of subsystems increases, the benefit of higher 
efficiency decreases. The resulting state of charge (SOC) pattern of the batteries 
is similar for the investigated systems, while the average SOC value is higher in 
highly modular systems. 
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1 Introduction 


Grid storage systems on megawatt scale play a vital role for the integration 
of renewable energies into electricity markets and grids. Several investigations 
focus on the development of optimized battery operation strategies [1-5]. For 
several reasons, existing grid storage systems usually consist of multiple batteries 
and inverters. For reasons of economies of scale, hardware manufacturers offer 
components in limited number of size classes, allowing lower production costs. 
Furthermore, the use of modular components supports the systems’ scalability. In 
addition, the size of modular systems can be changed over their lifetime. Moreover, 
modular systems avoid single points of failure, which leads to higher fault tolerance. 
Although the modularity of existing grid storage systems is well known, most of the 
analyses describe the storage system as a single battery combined with a single 
inverter [1-9]. Few studies are known that analyze the modularity of grid-connected 
battery systems and the related effect on system efficiency during operation and the 
influence of these energy losses on the operating strategy of the system [10, 11]. 

The consideration of this architectural aspect in the model-based analyses of 
battery operation provides a degree of freedom in optimizing the overall yield of 
a grid storage system. Figure | shows an example for a grid storage system with 
a high modularity. In this setup, three inverters and batteries are connected to one 
point of common coupling. 

Storage systems operated together with a renewable energy source are most likely 
not charged and discharged at their nominal power. Figure 2 shows the frequency 
of the operational inverter power over the course of a year for a large battery 
that is operated together with a wind farm with the purpose of supporting market 
integration of wind energy [14]. 

During 79% of the time, the system is in standby mode. During 5% of the time 
the system is operated at less than 50% rated power, while during 14% the power 
rating is between 80% and 100%. 

This suggests the importance of energy conversion efficiency not only in full- 
load operation, but also in standby mode and in partial-load operation. The inverter 
efficiency, however, is a nonlinear function in terms of power flow (Fig.3 with 


Fig. 1 Example for a grid storage system realized as a combination of three single battery inverter 
units (N=3) 
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Fig. 2 Frequency distribution of operational inverter power in our case study over one year 
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Fig. 3 System efficiency for different number of inverters (N) 


N= 1). In partial-load operation, the relative inverter losses are significantly higher 
in comparison to relative losses at full-load operation. Because of the possibility 
to use several inverter units, the efficiency is higher in high modularity systems. At 
low power flows, using one smaller inverter unit can limit the losses. Figure 3 shows 
the effect of number of inverters on the total efficiency curve. Here, we assume the 
power rating of the total storage system to be constant. Therefore, the power rating 
of the inverter and the constant losses scale with +- The impact on the efficiency is 
evident up to a number of inverters of 4. Constant losses become less relevant with 
the increase from 4 to 8 inverters. 

Figure 3 provides an overview of the influence on the efficiency of modular 
battery systems. For operating a system with higher modularity, with discrete 
battery-inverter pairs, the operating strategy is more complex. Without an optimiza- 
tion strategy, which combines knowledge about power transfer losses, the SOC of 
the batteries and requests from the grid, suboptimal operating states can be realized, 
which increases losses. Hence, the operating strategy has a significant influence on 
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the yield of a grid storage system. In this study we avoid choosing an individual 
optimization strategy by formulating an optimization problem. This optimization 
problem is formulated in terms of power flows and minimizes the power losses of 
the system. So, changes in the operating strategy, caused by increasing the number 
of inverter-battery pairs, are identified by solving the optimization problem. 

The target function optimizing the operating strategy of the system is purely 
technical driven, i.e. we optimize in terms of power losses. Since most of the 
business models for grid storage systems rely on the power in and outflows of the 
system for economic evaluation, an optimization in terms of losses is always of 
benefit for the system operator. 

As there is no direct connection between the batteries on the DC-side the energy 
management needs to take the state of charge of each battery into account. This 
increases the complexity of the energy management strategy. 

This study addresses the questions to which extent the segmentation of the system 
into multiple battery-inverter subsystems can reduce operating losses and aims 
at quantifying this effect. Therefore, we develop a mixed-integer linear program 
that represents the operational strategy of an energy management system of such 
a modular system as it could be applied in real applications. The mixed-integer 
linear program determines the system operation for meeting a given schedule for 
the whole system at the point of common coupling by minimizing inverter losses. 
As a result, the required energy to be fed into or stored from the grid is allocated 
over the different battery-inverter-subsystems. The resulting SOC of the batteries is 
investigated and ideas for further investigations are derived. 

Furthermore, the overall system topology has an influence on the power losses 
and the operating strategy. This study reduces the complexity by looking for 
identical battery-inverter pairs, which are coupled on the AC-side. A coupling on 
the DC-side has not been investigated in this study. 

This paper is structured as follows. Section 2 describes the mathematical model 
and presents the data chosen for the analysis. Section 3 presents the results and 
compares the efficiencies of systems with different degrees of modularity. Section 4 
draws conclusions for the design of grid-connected battery systems. 


2 Methodology 


2.1 Mathematical Model 
2.1.1 Power Flow Model 


The grid storage system is described as a set of nodes, representing the discrete 
time steps ¢ with duration of the time steps d, with transfer of energy from one 
node to another [12, 13]. In this model, each pair of inverter and battery represents a 
single storage node. The inverter is only described by its influence on the power 
transfer from the battery into the grid and vice versa. The SOC represents the 
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battery state of charge and interlinks the periods between each other, adding energy 
flows, both for energy charged and energy discharged during the previous period. N 
storage systems S; are connected to one point of common coupling. G’ represents 
an additional connection to the grid. It serves as a backup when power requests to 
the storage cannot be met and shall avoid infeasibility of the model. 

Hereinafter, a power flow from node A to node B is described as Ag and the 
efficiency of this transfer is described as na» (Ag). Per definition, transfer losses 
are assigned to the sink. 

The underlying power model is generic and initially independent of the technical 
implementation of the storage system. Different technologies might add different 
boundary conditions and parameters to the power flow. All parameters underlying 
the assumed technical realization are described in 2.1.3 to 2.1.5. 


2.1.2 Consideration of Losses 


The power transfer losses can be divided into constant, linear, and quadratic terms. 
This results in the common representation of the inverter efficiency (1): 


Pout = Pin — (Az + ba, - Pin + CAg' P>) (1) 


In case of battery storage systems, constant losses a4, have their origin in 
auxiliary power, e.g. for the battery management system, active cooling and thermal 
control of the battery cells, and other subsystems. These losses are independent 
from the status of the inverter. Some systems can reduce the constant losses during 
standby operation by switching off subsystems. We consequently assume that the 
systems considered in our analysis can be turned off, avoiding standby losses during 
operation, as it is shown to be possible in Munzke et al. [14]. 

Linear losses ba, are proportional to the rated power. They are mainly heat 
losses. Furthermore, we consider battery efficiency. This is reflected by losses that 
are proportional to the charging and discharging power. 

Quadratic losses c4, represent losses caused by nonlinear saturation effects in 
the inductance. In this work, the quadratic terms are not taken into account because 
non-linear losses cannot be clearly observed in all inverter systems [14]. 

Self-discharge of the battery cells, here represented as power loss of the storage 
path nss, is usually very low and thus neglected [14-16]. Moreover, we do not 
consider the power supply of the battery management system (BMS). 

We furthermore do not take into account losses which are independent of the 
degree of modularity, such as the power consumption of sensors or the climate 
control. Moreover, we do not consider transformer losses, neglecting the effect of 
additional power flows, introduced for ensuring model feasibility, on transformer 
losses. 

Inverter loss data is obtained based on a curve fitting approach of the efficiency 
curve of SMA’s 250 kVA inverter “Sunny Central” [17]. The battery loss parameters 
are approximated based on data measured in Munzke et al. [14]. 
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2.1.3 Parameters and Decision Variables 


Table 1 shows the exogenously given parameters used in the model. Table 2 lists the 
decision variables as well as their lower and upper bounds. Power values are always 
given in Watt. The penalty parameter p and two decision variables for additional, 
unplanned power flows to and from the grid (G7, and GL) are introduced for 
reducing deviations from the predetermined schedule and at the same time ensuring 
the model solvability. The same penalty parameter is applied both for charging from 


and discharging into the grid. 


Table 1 Description of parameters 


Parameter Description Unit 
Gy Maximum charge power W 
S ie Maximum discharge power W 
K Storage power capacity WwW 
N Number of invenen — D 
T Number of time steps per optimization = 

d h 

=0 = 

SOC; Battery state of charge at t = 0 Wh 
as;, Constant power losses for discharging (of ith storage) W 
AGs, Constant power losses for charging WwW 
bs, Linear inverter power losses for discharging % 
bes, Linear inverter power losses for charging % 
CSi Linear battery losses for discharging % 
CGs, Linear battery losses for charging % 
bs, Linear power losses due to self-discharge of ith storage % 
G! Power demand at point of common coupling W 
p Penalty for additional power flow = 


Table 2 Description of decision variables 


Gg Power flow grid to ith storage [0, 65] 

Sic Power flow ith storage to grid [0, Ss] 
SOC Battery state of charge of the ith storage [Wh] [0, « - £] 

YGs, Binary variable for charging of ith storage e {0; 1} 

Ys; Binary variable for discharging of ith storage e {0; 1} 

Gi Additional power flow to grid [0, Gy - N] 
Go Additional power flow from grid (0, Sic -N] 
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2.1.4 Target Function 


The target function minimizes all losses and penalizes additional power flows (2). 
The solution is obtained by carrying out one optimization for every hour, i.e. four 
15-min time steps at a time. 


T N 


min Y = > > (as VS + aGs, Yes, + (bs,, + Csia) Sig 
t=1 i=l (2) 


+ (bos, + cas) Gs, + (GL + Ge) +P) 


The target function determines the behaviour of the optimized power management 
strategy. In a modular system, the given power demand can be met by a subset of 
inverters. Given the example of two inverters in Fig.3, only one inverter is used 
if it can provide the requested power, thus limiting constant losses. When power 
demand exceeds the single inverter’s power rating, the second subsystem is used. In 
this case, there is no incentive to operate the modules at different power. Therefore, 
if charge and discharge power requests are sufficiently high, we expect a nearly 
equal distribution on the SOC. 


2.1.5 Constraints 


Two binary variables are introduced for charging and discharging (yG s; YSig = 
{0; 1}) in order to ensure that the storage system is not charged and discharged at 
the same time (3). 


YGs, + YSig S 1 (3) 
The power flow at the point of common coupling is defined by Eq. (4): 
G' =G; — Gb 
ee t t t (4) 
_ Sn (Gs, + (1 — bs, Sig — aSig Ysi) 
t=1 i=1 


Equation (5) describes the energy flow for each storage system. It is solved for 
each hour independently. Therefore, the SOC‘ is set as an external parameter. For 
the subsequent iterations, the SOC-value of the first time step is set as parameter 
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according to the solution of the optimization of the previous optimization. 


SOC! = (a — bs, )SOC!! + (1 — bas : cas, GS, 
(5) 


U = by, €5,,)Sig — Wig Vig Z 265, Vs, ) 4 


In addition, Eqs. (6) and (7) assure that battery charging or discharging is only 
realized if the binary variable is equal to 1. 


Gis, — OS Ys 0 (6) 
er, S 9 (7) 


Equations (8) and (9) set the boundaries for the power demand at the point 
of common coupling. These boundaries are calculated outside of and prior to the 
optimization. 


e 8) 


G' < GY. N (9) 


2.1.6 Implementation 


The model is implemented in MATLAB, using CPLEX as a solver. The solution is 
obtained on an hourly basis in 15-min resolution with a MIP gap of 0.5%. 


2.2 Data 


The battery schedule was calculated by Ried et al. [4]. It results from a planned 
operation of a 50 MW/100 MWh battery which is connected to a SOMW wind 
farm. The system participates in the German day-ahead and tertiary control reserve 
markets. The battery schedule is a time series in 15-min resolution and obtained by 
a mixed-integer linear program maximizing the contribution margin of the system. 
For this study, the schedule is scaled to a 2 MW/2 MWh battery system and used as 
the power demand at the point of common coupling G’. Since it does not account for 
detailed losses, the battery used for this study is scaled 30% larger. All assumptions 
are given in Table 3. 

Based on this data, the model is applied to six systems with varying degrees of 
modularity (Table 4). 
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Table 3 Parameter values 


Parameter Value 
Duration of one time step d 0.25h 
Number of time steps T 4 

Storage power capacity K 2.6 MW 
Battery SOC SOC!=? 1.3 MWh 
Maximum charge power Gy" 2.6 MW 
Maximum discharge power S; 2.6 MW 
Constant discharging losses as,_, 4.5% Gs 
Constant charging losses aG, 4.5% GS” 


Linear discharging losses (inverter) bs, 2.1% 


Linear charging losses (inverter) bg s; 2.1% 
Battery charging losses cGs, 2.5% 
Battery discharging losses cs, 2.5% 
Linear losses of self-discharge bs, ; 0% 
Penalty for additional power flow p 1015 


Table 4 Analyzed inverter numbers 


No. of inverters 


Battery size KWh 
an 2,600 
2 hw 1.300 

4 |6 650 
E 325 

16 |1625 162.5 
3200 [81.25 81.25 


3 Results and Discussion 


In this section, we compare the yearly losses and resulting battery operation for the 
different systems, and discuss potential implications for the system layout. 


3.1 Losses 


Figure 4 shows the losses of systems with varying degrees of modularity relative 
to the losses of the system with N=1. In accordance with the assumption shown 
in Fig.3, the losses scale with the number of inverter-battery subsystems. By 
increasing the number of inverters to 32, the operating losses can be reduced by 
38%. The gradual decrease of operating losses declines as the modularity increases. 
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Fig. 4 Yearly relative losses for different numbers of inverters 
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Fig. 5 Charge and discharge power in a single inverter system during operating mode 


This effect can be explained if the distribution of charge and discharge power 
over the course of one year for one and 32 inverters are compared. In a single 
inverter system, 50-70% of the charge and discharge power, neglecting standby 
mode, lies between 70% and 80% of the rated inverter power (Fig. 5). While 91% of 
the energy is charged at a rated power above 70%, 88% of the energy is discharged 
above 70% power rating. Power below 50% is requested during 29-31% of the time, 
corresponding to 5-6% of the energy flow. This is the mode of operation in which 
the individual inverter is less efficient and losses increase. 

As the number of inverters increases, the maximum power of each inverter 
decreases. Figure 6 shows the resulting distribution of charge and discharge power. 
As expected, the most likely power is the maximum power of the inverter, which 
corresponds to + of the maximum power of the single inverter setup. Only very few 
events occur, where low power is requested. 
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Fig. 6 Charge and discharge power in a system with 32 inverters during operating mode 
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Fig. 7 The distribution state of charge over one year for a realization with one (left) and 32 (right) 
battery systems 


3.2 Battery Operation 


A difference between the systems with one and 32 inverters is a different SOC-level 
during the course of the year. Figure 7 shows the distribution of the state of charge 
over the course of a year for systems consisting of one and 32 sub-systems. In the 
non-modular system, the battery is operated at an average 13% SOC, while in the 
high-modularity system the average SOC is 27%. While 86% of the SOC-values in 
the non-modular system range between 10-20%, 57% of the SOC-values are in this 
range in the system consisting of 32 inverters and batteries. Due to lower operating 
losses, less energy is wasted and higher SOC-values occur more often. 

Due to a larger number of batteries, there are relatively more downtimes of the 
batteries in high-modularity systems. While the system consisting of one inverter 
stands still during 79% of all periods, the sub-systems of the high-modularity 
systems are in standby-mode during 87% of the time. This suggests a higher 
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redundancy of the system consisting of 32 batteries, which could be beneficial in 
case of failure, especially when power output must be guaranteed. 


4 Conclusions 


In summary, we have performed a study on the operating losses of grid-connected 
battery storage systems on megawatt-scale consisting of several battery-inverter 
subsystems. Therefore, we applied a mixed-integer linear program determining the 
optimal operation of all sub-systems for a given schedule for the point of common 
coupling of a case study where a battery is operated together with a wind farm. 

The analysis shows that by increasing the system modularity, the frequency 
of operating points at low power decreases, resulting in a reduction of operating 
losses. The higher degree of freedom in modular systems thus allows following 
a given schedule at lower losses. By increasing the modularity from one to 
32 sub-systems, operating losses can be reduced by almost 40%. The effect of 
modularity on efficiency is most evident for systems consisting of few sub-systems. 
With increasing modularity, the gradual decrease of operating losses declines. The 
resulting state-of-charge of the batteries shows a similar operation pattern for the 
investigated systems. The avoidance of losses in high-modularity systems, however, 
leads to a shift towards higher SOC levels. 

We believe that this work motivates a modular layout of grid-connected battery 
systems. A higher efficiency should translate into lower operating cost. Moreover, 
the higher redundancy of modular battery systems can be advantageous particularly 
in applications that must ensure the system’s availability. The related investment 
cost is another aspect to consider when determining the system design. On the one 
hand, very large inverter and battery systems might translate into lower investment 
cost due to less balance of system components. A higher degree of modularity 
could thus be associated with higher costs. The opposite could also be the case, 
if the modularity enables the usage of standardized components produced in higher 
volumes, overcompensating the higher costs for the peripheral components. 

Further analyses could address the effect of considering monetary aspects within 
the target function. Moreover, the penalty parameter ensures model feasibility, 
but does not fully avoid additional power flows. Both the additional power flows 
and the choice of the penalty parameter could be further investigated. In order 
to penalize any deviation from the predetermined schedule, the same penalty 
parameter was applied both for charging from and for discharging into the grid. 
However, especially deviations leading to battery discharging might be less critical 
and sometimes even desirable from energy system perspective, possibly leading to 
positive revenues. This could be the case if the system participates in additional 
electricity markets. Consequently, the value of the penalty parameter could be 
different for positive and negative deviations and should be subject to further 
research. 
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Maybe the most interesting potential for further research lies in the operation of 
the different batteries, the effect on battery ageing and implications on technology 
choice and system design. Our results show higher average SOC-levels of the 
batteries in high-modularity systems. The average SOC has an impact on the 
calendar ageing, so that higher SOC-levels could lead to reduced battery lifetime, at 
least for most technologies. Another aspect worth investigating is the effect of higher 
C-rates on battery degradation. Because of allocating the same requested power to 
fewer batteries with lower capacities, the C-rates might be higher in high-modularity 
systems. The modular system architecture might thus be more suitable for some 
battery technologies than for others. The selection of one or more appropriate 
battery technologies could positively influence battery lifetime. 

Lastly, the battery management system determining the operational strategy 
could be adapted in order to take quadratic losses into account or vary operating 
patterns, SOC levels, and power requests for different storages. Moreover, asym- 
metric system topologies would allow for stressing the batteries at different depths 
of discharge and C-rates, even when running the same schedule. This might be an 
argument for hybrid storage systems consisting of different battery technologies 
with different characteristics and furthermore improve cycle ageing of the batteries, 
another important factor for battery degradation. 
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